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SECTION  I 
INTRODUCTION 

The  following  investigation  deals  with  constitutive  relations 
employing  generalized  derivatives  that  relate  stress  and  strain  in 
viscoelastic  materials,  and  the  solution  techniques  for  the  resulting 
equations  of  motion  for  structures  incorporating  viscoelastic  components 
to  damp  vibratory  motion.  The  use  of  generalized  derivatives  of 
fractional  order  in  stress-strain  constitutive  relations,  first  suggested 
by  Caputo  (Reference  1),  may  be  viewed  as  an  extension  of  the  standard 
model  for  a  linear,  viscoelastic  material. 


The  standard  viscoelastic  model  for  a  uniaxial  constitutive  relation 
is  (Reference  2) 


K 


o  (t) 


,k 

v  i  d  o 

l  bi-  ~ 


k=  1 


V(t> 


J 

l 

j  =  l 


E. 

J 


c 

df1 


(U 


The  viscoelastic  constitutive  relation  employing  generalized  derivatives 
of  fractional  order  will  be  taken  to  be 


o(t) 


K  t 

l  bkU 

k  =  l  k 


;[o(t)J 


E0c(t) 


+  I  Fi.Daj[t(t)]  (2) 

j  =  l  -1 


where  the  generalized  derivative  operator  of  real  order  a  is  defined  by 


n“[x(t)]  = 


1  d  f 

m-a)  cTt  • 


X  (t) 
(t-Tla 


dt 


0  <  a  <  1 


(3) 


The  generalized  derivative  constitutive  relation  (Equation  2)  may  be 
viewed  as  an  extension  of  the  standard  model  (Equation  1)  in  the  sense 
that  the  derivatives  are  no  longer  limited  to  being  of  integer  order. 
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The  use  of  this  generalized  derivative  constitutive  relation  in 
modeling  the  response  of  viscoelastic  materials  will  be  seen  to  have 
several  advantages  over  present  methods.  The  major  drawback  of  the 
standard  viscoelastic  model  (Equation  1)  is  that  a  large  number  of  terms 
are  often  required  to  describe  a  material  adequately.  The  use  of 
derivatives  of  other  than  integer  order  in  the  constitutive  relation  will 
be  seen  to  produce  satisfactory  models  with  very  few  parameters.* 

Because  of  the  large  number  of  terms  required,  the  standard  model 
(Equation  1)  often  becomes  too  cumbersome  to  manipulate.  Consequently, 
an  alternative  known  as  the  "complex  modulus  method"  has  been  developed. 
In  the  complex  modulus  method,  measured  values  of  E*(w)  (Equation  4)  are 
used  as  a  discrete  approximation  of  the  function  E*(oj).  In  the  transform 
domain,  the  general  viscoelastic  constitutive  relation  is 

a*M  =  E*(u)e*(u)  (4) 

E*(w)  is  measured  for  different  frequencies  of  motion  (Reference  4),  uk  , 
which  produces  a  set  of  discrete  values  of  the  modulus,  E*(u).j ),  over  the 
frequency  range  of  interest.  These  discrete  values  of  E*(w)  are  sub¬ 
stituted  into  the  transformed  equations  of  motion  of  a  viscoelastic 
material  to  produce  values  of  the  transform  of  the  response  at  discrete 
frequencies.  The  inverse  transform  of  the  response  is  evaluated 
numerically  to  produce  the  time  history.  The  major  drawback  of  this 
method  is  the  arduous  task  of  calculating  the  inverse  transform  for  every 
point  in  time  at  which  the  value  of  the  response  is  required.  The  use  of 
the  generalized  derivative  constitutive  relation  will  do  a way  with  the 
need  for  numerical  approximations  in  the  frequency  domain. 

An  elementary  form  of  the  "complex  modulus"  method,  obtained  by 
representing  the  transform  of  the  modulus  by 

E*(«0  =  Eq(1  +  insgn(w))  (5) 


*In  Section  V,  a  generalized  derivative  constitutive  relation  for  the 
elastomer  3M-467  is  presented.  The  relation  characterizes  the  material's 
properties  over  four  decades  of  frequency  with  three  parameters. 
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u  >  0 

a)  =  0  (6) 

u  <  0 

and  often  known  as  "structural  damping"  (Reference  5),  is  valid  only 
for  sinusoidal  stress  and  strain  in  the  material.  Crandall  has  shown 
that  the  response  of  a  harmonic  oscillator  with  "structural  damping" 
for  impulsive  loading  is  non-causal  (Reference  6);  that  is,  the  time 
response  of  the  oscillator  occurs  before  the  loading. 

Milne  (Reference  7)  has  proposed  several  modifications  of  the 
imaginary  part  of  the  modulus  given  in  Equation  7  to  produce  a  causal 
response.  Unfortunately,  neither  the  modified  modulus  nor  the  one  given 
in  Equation  7  is  particularly  suitable  for  transient  (broad-band)  response 
of  viscoelastic  materials,  because  they  do  not  account  for  the  frequency- 
dependent  stiffness  typically  encountered.  Caputo  (Reference  8)  observed 
that  a  single  term,  generalized  derivative,  constitutive  relation  of 
the  form 

aft)  =  E1Dai[E(tj]  0  <  cxl  <  1  (7) 

produces  frequency-dependent  stiffness  and  damping  and  a  loss  factor*, 
n,  that  is  frequency  independent. 

a,  it 

n  =  tan  (8) 

it 

Caputo' s  work  with  generalized  derivative  constitutive  relations 
focuses  primarily  on  the  propagation  of  waves  in  geological  formations. 

One  of  Caputo's  earliest  papers  (Reference  9)  was  on  generalized 
independent  loss  factors  or  equivalently  frequency- independent  resonant 
qualities,  Q.  This  work  was  followed  by  a  book  (Reference  10)  in  which 
Caputo  dealt  with  the  propagation  of  impulsive  plane  waves  and  the 


*The  loss  factor  in  a  linear  material  is  the  ratio  of  imaginary  to  the 
real  part  of  the  modulus. 


j  1, 

sfin(w)  =  0, 

(-1  , 
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free  vibration  of  spherical  strata  using  generalized  derivatives  in  the 
constitutive  relations  of  the  media.  In  the  last  chapter  of  the  book 
and  in  later  papers  with  Minardi,  Caputo  compares  the  generalized 
derivative  relations  with  experimental  observations  of  the  properties 
of  some  media:  "some  metals,  glasses  and  the  earth."  (References  11,  12) 
In  1974  Caputo  proposed  a  generalized  derivatives  viscoelastic  con¬ 
stitutive  relation  of  the  form  (Reference  13) 

o(t)  =  nDa+n[c(t)]  ,  0  <  a  < 1  ;  n  =  0,1,2,...  (9) 

where  the  generalized  derivative  operator  was  defined  by 

Da+n[x ( t } ]  s  /l  ~~T  dt  (,0> 

v  0  dt  (t-r) 

This  definition  differs  slightly  from  the  generalized  derivative  used 
in  this  investigation  (Equation  3).  Caputo  used  the  relation  (Equation  9) 
to  determine  the  response  of  a  uniformly  driven  infinite  viscoelastic 
layer  and  investigated  the  hysteresis  behavior  of  the  constitutive 
relation  (References  14,  15).  Recently,  Caputo  suggested,  but  gave  no 
application  for,  a  constitutive  relation  of  the  form  (Reference  16) 

DP[a(t)]  =  nDa  [  e  (t )  ]  (11) 

which  is  the  harbinger  of  the  general  constitutive  relation  (Equation  2) 
which  is  to  be  used  in  this  study. 

All  of  Caputo' s  work  rests  on  a  continuum  formulation  of  the  equations 
of  motion.  Unfortunately,  a  continuum  formulation  of  the  equations  of 
motion  for  many  structures  of  engineering  interest  is  not  practical.  As 
a  result,  a  discrete  formulation  of  the  equations  of  motion  is  adopted, 
based  on  assumed  displacement,  finite-element  methods.  A  solution 
technique  developed  for  the  motion  of  structures  incorporating  visco¬ 
elastic  materials  modeled  with  generalized  derivatives  is  developed  as  an 
extension  of  the  solution  technique  developed  by  Foss  for  non-proportional 
viscous  damping  (Reference  17). 
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In  the  sections  to  follow,  constitutive  relations  using  generalized 
derivatives  are  developed,  and  their  applications  and  limitations  are 
considered.  Some  experimental  results  are  presented  and  found  to  show 
that  a  constitutive  relation,  with  parameters  determined  from  the 
response  to  sinusoidal  loading,  predicts  very  well  the  response  of  one 
typical  elastomeric  material  to  an  impact  loading.  Finally,  the 
formulation  of  multi -degree  of  freedom  systems,  necessary  for  large- 
scale  structural  analysis,  is  considered.  Special  solution  methods, 
necessary  for  such  applications,  are  developed. 
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SECTION  II 

A  BRIEF  OVERVIEW  OF  GENERALIZED  DERIVATIVES 

Before  construction  of  generalized  derivative  constitutive  relations, 
it  is  appropriate  to  introduce  the  properties  of  generalized  derivatives 
relevant  to  the  following  investigations.  Of  particular  interest  are  the 
form  of  the  Laplace  and  Fourier  transforms  of  generalized  derivatives, 
and  the  results  of  repeated  differentiation  of  fractional  order. 

First  and  foremost,  the  generalized  derivative  is  a  linear  operator. 
Da[Xl(0  +  x2(t)]  =  D“[x1Ct)]  +  Da[x2(t)]  (12) 

This  property  follows  directly  from  the  definition  (Equation  3) 

2  <3) 

To  put  the  definition  into  a  form  in  which  the  calculation  of  its 
Laplace  transform  is  straightforward,  one  first  performs  a  change  of 
variable 

t  =  t-n  (13) 


which  results  in 


'  T7lW  3F  l  dn 


(14) 


Using  Leibnitz's  rule  to  differentiate  the  integral  produces 

"“[’‘(03  ’  TTT^T  4  4  *  rTTW  <,5> 


o  n 
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After  taking  the  Laplace  transform  of  Equation  14,  the  transform  of  the 
generalized  derivative  of  order  a  of  the  function  x(t)  is  seen  to  be 

L[Da  [x  (t )  j  3  =  -J—  •  (sL[x(t)]  -  x  (0) )  +  (16) 

s  s 


which  simplifies  to 

L[D"[x(tj]]  =  saL[x(t)3  (17) 

where 

00 

L[x (t )  ]  =  /  x(t)e"st  dt  (18) 

o 

Notice  that  the  Laplace  transform  of  a  generalized  derivative  of  order  a 
of  a  function  is  equal  to  sa  times  the  transform  of  the  function. 

Under  certain  conditions  a  similar  property  of  generalized  derivatives 
is  true  for  Fourier  transforms. 

F[Da[x(t)]]  -  (i„)a  F[x  (t )  ]  (19) 

where 

F[x  (t )  ]  s  /  x(t)e'iu,t  dt  (20) 

-  ao 

The  conditions  are,  first,  that 

x(t)  =  0  for  t  <  0  (21) 

in  which  case  the  Fourier  transform  becomes 

F[x (t ) ]  =  /  x(t)e'iwt  dt  (22) 

o 


and,  second,  that  the  integral  in  Equation  22  exists.  Note  the  parallel 
form  of  Equations  17  and  19.  Both  relations  were  used  by  Caputo 
(Reference  18). 
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A  useful  property  of  generalized  derivatives  is  that  the  generalized 

derivative  of  order  ,  of  the  generalized  derivatives  of  order  of  a 

function  is  the  generalized  derivative  of  order  a-]  +  ct^  of  the  function. 

Again  using  operator  notation,  the  property  is 

a.  a9  ai+a? 

D  1 [D  2[x(t)]]  =  D  1  2[x(t)]  (23) 


Notice  that  the  definition  of  generalized  differentiation  given  in 
Equation  3  is  restricted  to  fractional  order  a  less  than  one.  If  a  is 
one  or  greater  in  Equation  3,  the  integral  contains  a  non-integrable 
singularity*.  The  definition  of  a  generalized  derivative  of  order  3, 
where  6,  where  6  >  1  and  6  =  m  +  a  where  m  is  the  largest  integer  not 
exceeding  3,  is 


D"+“  [x(t)] 


rTT^T 


m+1 


dt 


m+1 


t 

/ 


o 


x(t) 

(t-T) 


a 


dt 


(24) 


*0n  the  other  hand,  if  a  is  zero  in  Equation  3,  the  relation  is  clearly 
valid  and  follows  from  the  fundamental  theorem  of  calculus. 
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SECTION  III 

THE  BASIC  GENERALIZED  CONSTITUTIVE  RELATION 


In  this  section,  the  basic  generalized  derivative,  viscoelastic 
constitutive  relation  is  presented  and  some  aspects  of  its  behavior  are 
established.  The  basic  generalized  derivative  constitutive  relation  is 


or 


o(t)  -  XDa[  e  (t )  ]  0  <  a  <  1  (25) 


oft) 


X  d  fl  e  (n) 

l 


dn 


(26) 


Since  the  generalized  derivative  is  a  linear  operator,  this  relation 
is  suitable  only  for  the  linear  approximation  of  a  material's  properties. 
This  linear  constitutive  relation  satisfies  many  of  the  presently  accepted 
constraints  on  viscoelastic  constitutive  relations. 


In  particular,  it  represents  a  material  with  fading  memory.  To 
demonstrate  this  claim,  it  is  necessary  to  put  the  constitutive  relation 
in  a  different  form  using  a  change  in  variable. 

n  =  t-T  (27) 

and  again  using  Leibnitz's  rule  to  differentiate  the  resulting  integral 
produces 


oft) 


I'D -a) 


I1  4 


a 

at 


eft-t) 


dr 


x  E  (0  ) 
r(l-a)ta 


(28) 


or 


t 

aft)  =  /  G(T)e(t-T)dx  +  Gft)e(O) 

o 


(29) 


where 


Gft) 


X 

rfl-a)ta 


(30) 


1 
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The  constitutive  relation,  (Equation  29),  represents  a  viscoelastic 
material  with  a  fading  memory.  A  material  is  said  to  have  a  fading 
memory  if  its  relaxation  modulus,  G(t),  goes  to  zero  monotonically  as  t 
increases  (Reference  20).  Notice  that  G(t)  in  Equation  30  does  in  fact 
go  to  zero  monotonically. 

Since  the  material  has  a  memory,  the  value  of  the  stress  at  time  t 
is  dependent  on  the  entire  strain  history  until  time  t.  To  ensure  that 
the  constitutive  relation  produces  a  stress  that  is  dependent  on  the 
entire  strain  history  until  t,  time  zero  must  be  chosen  before  or  at  the 
onset  of  the  initial  strain. 

Consequently, 

e(t)  =  0  for  t  <  0  (31) 

and  the  only  way  that 

e (0 )  ^  0  (32) 

is  if  the  strain  history  is  discontinuous  at  t  =  0.  According  to 
Gurtin  and  Sternberg  (Reference  21),  the  constitutive  relation  as  shown 
in  Equation  29  is  in  the  correct  form  to  handle  discontinuous  strain 
histories.  Notice  that  a  step  discontinuity  in  the  strain  history  at 
t  =  0  produces  a  stress  history*  that  is  singular  at  t  =  0. 

When  the  strain  history  is  a  continuous  function  of  time,  and  zero 
for  negative  times,  the  constitutive  relation  given  in  Equation  29  reduces 
to 

t  X 

o(t)  =  /  G(T)e(t-x)dT  ,  GCO  =  - r  (33) 

o  r(i-a)ta 

This  relation  satisfies  three  out  of  four  of  Pipkin's  restrictions  on 
viscoelastic  constitutive  relations  (Reference  22).  The  first  restriction, 
that  the  stress  be  an  odd  functional  of  strain  rate,  is  satisfied. 

♦The  Voigt  viscoelastic  model  displays  this  same  property. 
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The  second  restriction,  that  G(t)  go  to  zero  as  time  increases,  is  also 
satisfied,  which  is  in  keeping  with  the  fading  memory  property.  The  third 
restriction,  also  satisfied,  is  that  the  kernel,  the  relaxation  modulus 
G(t),  be  a  function  and  not  a  distribution.  Pipkin's  fourth  restriction, 
not  satisfied  by  the  generalized  derivative  constitutive  relation,  is 
that  G(t)  be  of  negative  exponential  order. 
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SECTION  IV 

GENERALIZED  DERIVATIVE  CONSTITUTIVE 
RELATIONS  FOR  VISCOELASTIC  MATERIALS 

The  task  at  hand  is  to  use  the  basic  generalized  derivative  con¬ 
stitutive  relation  just  presented  as  the  building  block  for  constitutive 
relations  that  model  the  frequency-dependent  moduli  of  viscoelastic 
materials.  The  moduli  of  viscoelastic  materials  are  complex  numbers 
where  the  real  anc  imaginary  parts  are  functions  of  the  frequency  of 
motion. 

A*(w)  =  A'(u)  +  iX"(w)  (34) 


U*(w)  =  u'(u)  +  i  p  "  ( a) ) 


(35) 


The  moduli  are  defined  as  the  transforms  that  relate  the  transforms  of 
stress  and  strain 

0*n(w)  =  6mnX* (w)e* (u)  +  2y*  (u) c*n (u)  (36) 

where  <5  is  the  Kronecker  delta  and  e*(w)  is  the  transform  of  the 
mn 

dilatation  strain 

e*(w)  =  cii  (“)  +  e22^  +  e33^w-*  (37) 


A  useful  property  of  the  moduli  is  that  their  values  at  frequency  w0 
relate  sinusoidal  stress  and  strain  of  frequency  coQ  in  the  material. 


°mn(t)  =  6mnX*(a,o)e0expCiu,ot]  +  ^oUmn  expfi^t]  (38) 

o 


Consequently,  one  can  measure  the  values  of  u*(u>)  and  A*(w)  at  discrete 
frequencies  of  sinusoidal  motion.  As  a  result,  the  frequency  dependence 
of  the  moduli  can  be  determined  experimentally. 


Typically,  a  viscoelastic  material  at  constant,  uniform  temperature 
has  moduli  that  vary  with  the  frequency  of  motion  as  indicated  in 
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Figure  1  (Reference  23).  At  low  frequencies  (the  rubbery  region)  the 
real  part  of  the  modulus  is  relatively  constant,  while  the  imaginary 
part  of  the  modulus  increases  with  increasing  frequency.  At  intermediate 
frequencies  (the  transition  region)  both  the  real  and  imaginary  parts  of 
the  modulus  increase  with  increasing  frequency  and  the  rate  of  increase 
of  the  real  part  slowly  overtakes  the  rate  of  increase  of  the  imaginary 
part.  At  high  frequencies  (the  glassy  region)  the  imaginary  part  of  the 
modulus  decreases  with  increasing  frequency,  and  the  real  part  of  the 
modulus  is  relatively  constant. 

The  generalized  derivative  constitutive  relations  presented  here  are 
of  two  types.  The  first  type  are  those  relations  intended  to  model  the 
viscoelastic  behavior  of  the  material  in  the  rubbery  and  transition 
regions.  For  brevity,  this  type  of  model  is  referred  to  as  the  RT  model. 
The  second  type  are  those  relations  intended  to  model  the  behavior  of 
the  material  in  the  rubbery,  transition  and  glassy  regions.  This  type  of 
relation  is  referred  to  as  the  RTG  model.  Since  the  RTG  model  may  be 
viewed  as  a  generalization  of  the  RT  model,  the  RT  model  is  considered 
first. 

The  RT  model  for  an  isotropic,  homogeneous,  linear  viscoelastic 
material  is 


o  „(t) 
mn 


6  (X 
mn  v  o 


J  a. 

+  l  X  D  J)e(t) 
3  =  1  3 


where 


(39) 


0  <  Qj  <  1 


(40) 
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0  "  “t  "  1 


(41) 


e(t)  =  CjjU)  +  E22^t'*  +  c33^t^  (42) 


Note  that  the  RT  model  portrays  the  stress  as  a  combination  of  elastic 
stresses,  proportional  to  the  positive,  real  parameters  XQ  and  and 
viscoelastic  stresses,  proportional  to  the  positive,  real  parameters  A. 
a  nd  u , . 

*v 

To  establish  the  frequency  dependence  of  the  moduli  in  the  RT  model, 
one  takes  the  Fourier  transform  of  the  constitutive  relation. 


V„<“>  ■  Smn(lo  ♦T/jO"')  3>  '*<“> 


2  U°  *Ji"«  <i":i  ^  e"n(“:i 


(43) 


Expressing  the  moduli  in  terms  of  their  real  and  imaginary  parts  produces 


*  r  >. 
A  (to) 


a  .  IT  x  • 

3  p  nc _ 3 


(A-  +  y  A  .  to  J  eos-7-J— ) 
0  j=l] 


J  a .  Ha • 

+  i  l  A  .  w  ^  S  in— 

j“l  3  Z 


(44) 


L  a  IIo . 

V  (u)  *  (un  +  1  W.w  COS  -y— ) 

°  £=1  1  *• 


.  y  ai  .  n% 

1  l  u  U)  sm 
1  =  1 


(45) 
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Observe  that  for  low  frequencies  of  motion,  defined  by 

<<  1 


,  J  a . 

-  I  J 


‘o  jii  1 


and 


—  I 

vo  i=l 


<<  1 


(46) 


(47) 


the  moduli  in  the  RT  model,  (Equations  44,  45),  have  essentially  a 
constant  real  part  and  an  imaginary  part  that  increases  with  increasing 
frequency,  similar  to  the  properties  of  a  viscoelastic  material  observed 
in  the  rubbery  region  (Figure  1). 


At  intermediate  frequencies  of  motion,  defined  by 

,  J  a . 

—  I  J  1-1 

o  j  =  l  J 


and 


1 


p 


:  l 


(48) 


(49) 


the  real  and  imaginary  parts  of  the  moduli  in  the  RT  model  are  increasing 
with  frequency,  similar  to  the  properties  of  a  viscoelastic  material 
observed  in  the  transition  region. 


It  is  evident  that  the  RT  model  does  not  properly  account  for  the 
properties  of  a  viscoelastic  material  at  high  frequencies,  defined  by 


1  r  ,  cx. 

T-  l  1;“  J 
0  j  =  l  J 


>>  1 


(50) 
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and 


1 


,  td  “  !. 


>>  1 


(51) 


The  real  and  imaginary  parts  of  the  modulus  are  predicted  by  the  IT  model 
to  increase  indefinitely  with  frequency.  In  the  glassy  region,  however, 
the  real  part  of  the  modulus  of  viscoelastic  materials  is  typically 
constant  (Figure  1),  and  the  imaginary  part  is  decreasing  with  increasing 
frequency.  This  discrepancy  motivates  the  construction  of  a  more  complex 
model  which  accounts  for  material  properties  in  the  glassy  region,  the 
RTG  model . 


The  RTG  model  for  an  isotropic,  homogeneous,  linear,  viscoelastic 
material  is  defined  to  be 


K  6,  P  IT 

(1  +  l  a,  D  k)  (1  ♦  l  b  D  P) 
k=l  k  p«l  p 


mn 


(t) 


mn 


(1 


p  e 

•  [  b  C  p) 

S  p  ’ 
p=l  v 


J  a- 

+  l  X.D  J)  e(t) 
3-1  J 


k  e,  l 

2  (1  *  l  a,D  k)  (v  +  l 
k=l  K  °  i  =  l 


p,d  *) 


(52) 


Again,  the  frequency  dependence  of  the  moduli  in  the  model  is  observed 
by  taking  the  Fourier  transform  of  the  constitutive  relation.  After 
some  algebraic  manipulation  of  the  transformed  relation,  the  result  is 


°mn 


M 


(X  ♦  l  X  (iu>)aj) 

5 _ _ jp* - s -  e  (w) 


mn 


Si 


(1+1  a.(iu)  ) 
k=l  K 
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where 


- p - : — 

(1  +  I  bn(iu)BP) 

P-1  P 


X  (ui) 


(x  +  l  x  Cico)  J) 
0  3-1  3 

P  §7“ 

(l  +  [  b  (iu)  k) 

P-1  P 


u*(«) 


Cw  +  l  v, (iw)  ; 

°  t=i  s- 

- p - 1 — 

C  1  +  I  b  (id,)  P) 

P-1  P 


0  <  aj  ’°Vek,ep  <  1 


(53) 


(54) 


(55) 


(56) 


and 


X„ ,  X  . 
0  J 


uo>ut*ak> 


If  the  parameters  a^ 


and  b 

P 


are  chosen  to  be  small,  such  that 


(57) 


and 


X* 

II 

ak  C^w) 

<  < 

>• 

o 

+ 

M  C — -3  CL_ 

h-» 

x j  (id,)  i 

P 

l 

P-1 

6 

bp  (id,)  P 

<  < 

L 

Wo  +  * 

l=l 

Wj  (iw) 

(58) 


(59) 
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for  frequencies  in  the  rubbery  and  transition  regions  of  the  material, 
the  RTG  model  behaves  like  the  RT  model  in  the  rubbery  and  transition 
regions  of  the  material. 


It  is,  however,  the  presence  of  the  terms  a.  and  b  which  enables 

K  P 

the  RTG  model  to  account  for  the  properties  of  viscoelastic  materials 
at  high  frequencies;  i.e.,  in  the  glassy  region.  Note  that,  if  the 
largest  values  of  and  are  the  same,  and  if  the  largest  values 
and  6p  are  the  same,  then  at  high  frequencies  A*(w)  and  p*{w)  have  real 
parts  that  become  constant  and  imaginary  parts  that  decrease  with 
increasing  frequencies,  as  is  characteristic  of  a  viscoelastic  material 
in  the  glassy  region. 


An  important  property  of  the  RT  and  RTG  models  is  that  they  satisfy 
the  "elastic-viscoelastic  correspondence  principle."  (Reference  4) 

The  correspondence  principle  states  that  the  Laplace  transform  of  the 
stress  response  of  a  viscoelastic  material  can  be  constructed  from  the 
Laplace  transform  of  the  response  of  an  elastic  material  by  replacing 
the  elastic  constants,  A  and  u,  in  the  elastic  response  by  the  Laplace 
transforms  of  viscoelastic  moduli,  A*(s)  and  u*(s).  The  principle  holds 
when  the  transform  of  the  elastic  stress-strain  constitutive  relation 


* 

a 

mn 


(s) 


6  Ae 
mn 


(5) 


■  v  e 


mn 


(s) 


(60) 


can  be  used  to  construct  the  transform  of  the  viscoelastic  stress-strain 
constitutive  relation  by  replacing  the  elastic  constants  with  the 
transforms  of  the  viscoelastic  moduli.  Thus,  the  viscoelastic  constitutive 
relation  must  be  of  the  form 


o  (s)  =  6  A*  (s )e*  (s) 

mnv  ;  mn  v 


2w‘(s)e;n(s) 


(61) 
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The  Laplace  transforms  of  the  RT  and  the  RTG  models  are  of  the 
general  form  given  in  Fquation  61.  The  Laplace  transform  of  the  RTG 
model  is 


* 

°mn 


(s) 


Smn  - h - F~ 

(  1  +  l  a.s  k) 
k=l  K 


P 

(  i  +  I  V  P) 

P=1  p 


e*(s) 


(62) 


and  the  Laplace  transform  of  the  RT  model  is 


Another  important  property  of  the  RT  and  RTG  models  is  that  they  can 
be  constructed  to  be  causal  in  the  sense  that  the  response  (stress)  does 
not  occur  before  the  input  (strain).  The  stress  response  is  zero  for 
negative  time  if  its  Laplace  transform  is  analytic  in  the  right  half  s 
plane.  This  condition  on  the  transform  of  the  stress  is  met,  for  the 
class  of  strain  histories  having  transforms  that  are  analytic  in  the 
right  half  s  plane,  when  the  branch  cuts  of  s°^,  s^,  s^k  and  sBP  are 
along  the  negative,  real  s  axis,  and  y*(s)  and  A*(s)  have  no  zeros  in 
the  right  half  s  plane.  Since  the  stress  is  zero  for  negative  time, 
the  stress  cannot  anticipate  the  strain  that  begins  at  time  zero. 

Laplace  transforms  are  also  useful  in  determining  the  hysteresis 
predicted  by  the  RT  and  RTG  models.  For  a  sinusoidal  strain  history 
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of  frequency  u>0  starting  at  time  zero, 
is 


the  transform  of  the  stress  history 


e,nnCt) 


Emn  “"V 
o 


(64) 


* 

o 

mn 


Cs) 


(6  >. 

mn 


Cs)e0  *  2u*(s)e„„  ) 


mn 


S  2  +  ^2 


(65) 


Evaluating  the  inverse  transform  in  the  same  manner  as  in  Section  X,  the 
stress  time  history  is  found  to  be 


o 

mn 


(t) 


i 

71 


( 6 

mn 


(ito  )  e 
o  o 


(i 


w^  )  e 

o  mn 


lu>  t 

e  o 


iw  t 


0  6°  *  2R* (*iw0) Gmn  )  e  "o 


+  Im 


v  f 


e0  *  2u*(re'llr)e  , 

o 


we 

o 


-rt 


u2+r2 


dr 


(66) 


As  time  increases,  the  integral  term  in  Equation  66  goes  to  zero  and 
the  sinusoidal  terms  dominate  the  stress  time  history.  So,  for  a 
sinusoidal  strain  history,  the  stress  eventually  becomes  sinusoidal 


as  wel 1 . 


Using  Euler's  formula,  numerous  trigonometric  manipulations,  and 
the  observation  that  the  two  sinusoidal  terms  are  conjugates,  the  stress 
for  large  time  may  be  evaluated  as  components  in-phase  with  the  strain, 
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proportional  to  sina^t,  and  components  out-of-phase,  proportional  to 
coswQt.  The  resulting  expression  for  the  stress  using  the  RTG  model  is 


aft)  = 
mn 

t>>0 


6  e 
mil 


J  a 


:  A.  (lO  )  X  +  V  X.u  i  "“j 
o  1  o  o  j  o  Jcos  — 2*- 


J  K 


K  &  .  J  K  a.+&, 

+  Xo  JjVo  cos  +  Vo  3  cos-f  (aj'&k) 


2  e_  A0(w  ) 
mn  2  o 
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L  a.  ira.  P  B  irB 
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o  i  o  IT  ojti  p  o  t 

1=1  p=l  v 


L  P  a  +B 
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+  y  y  p  b  (o 

1=1  p=i  p  p  ° 


sinu^t 
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K  S,  ttB-k  J  K  a  • +Bv. 

+  sinT +  jii  jJiXjV°  sinT(aj'ek) 
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L  P  a  +B_  _  -  - 

Jj  sinT  CVBp> 


COSw  t 

o 


(67) 


22 


AFML -TR-79-41 03 


where 


Al(uo) 


1  +  l  ak‘^iuo^ 
k=l  K  ° 


a2(uo) 


P  ** 

1  +  l  b  •  (iw  )  P 
p=l  P  0 


A  more  compact  form  of  the  expression  for  the  stress  is 


o  (t) 
mn 

t  >  >0 


6  F1 (u  )e  sinu  t  +  F-(w  )e  sinw  t 
mnloo  o  2^o  mnQ  o 


+  6  F,  (uj  )e  cosu  t  +  F.  (u>  )e  cosid  t  (70) 

mn3vo  o  o  oJ  mn  o 

o 

Under  conditions  of  uniaxial  stress  and  strain.  Equations  64  and  70  may 
be  recognized  as  the  parametric  equations  of  an  ellipse;  thus,  it  is 
clear  that  the  RTG  model  predicts  the  existence  of  a  hysteresis  loop  and 
the  loop  is  el  1  iptical . 


The  loss  factor  associated  with  the  hysteresis  loop,  the  ratio  of 
the  energy  dissipated  during  a  cycle,  D,  to  the  peak  strain  energy 
stored  during  a  cycle,  Umax,  is  a  parameter  often  used  to  characterize 
the  ability  of  viscoelastic  materials  to  damp  vibratory  motion. 


The  energy  dissipated  per  cycle  is 


2  v 

3  3  f  £“ 

D  =  l  f  o  (t)c  (t) 

v.-i  I  mn  mnv 

n=l  m=J  J 
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where  —  is  the  period  of  the  motion.  Using  Equations  64  and  70  to 
o 

evaluate  dissipated  energy  yields 


3  3 

d  =  7i  f  ,  ( w  )  e 2  +  y  y  F  .  ( w  )  c 

3^o  o  L.  4  o 

n=l  m=l 


2 

mn 


(73) 


The  peak  energy  stored  during  a  cycle  is 


U 


3  3 

max  =  11 

n=l  m=l 


Emn  Emn  ^ 

a  dc 
mn  mn 


E  =0 
mn 


2  a* 


o 

3  e  c 
mn  mn 


where  am„  are  the  stresses  in-phase  with  their  respective  strains,  c 
mn  r  mn 

o 

sinwQt.  Again  using  Equations  65  and  66,  the  peak  energy  stored  during 
a  cycle  is 


max 


1 

T 


Fl(“o)eo  * 


[  I  f,(»„)e2 


L  1  2 

n=l  m=l 


o'  mn. 


(75) 


The  resulting  loss  factor  is  seen  to  be 


n  = 


F3(uJo)eo  +  £  ^F4(“o)e 

n= 1  m=l 


2 

mn. 


Fl(%)eo  +  J.  J  F 2(wo)emn 
n=l  m=l  t 


(76) 


Notice  that  the  form  of  the  expressions  for  the  hysteresis  behavior 
and  loss  factor  of  the  RT  model  is  identical  to  that  of  the  RTG  model. 
This  follows  from  the  observation  that  the  RT  model  is  a  special  case  of 
the  RTG  model  for  =  0,  k  =  1,  2,  ...,  K,  and  bp  =  0,  p  =  1,  2, 

P. 
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In  summary,  RT  and  RTG  models  satisfy  the  elastic  viscoelastic 
correspondence  principle.  Conditions  necessary  to  ensure  that  the  stress 
does  not  anticipate  the  strain  have  been  developed.  In  addition,  both 
models  predict  the  existence  of  stress-strain  hysteresis  effects  and  the 
resulting  hysteresis  loops  are  elliptical.  Most  significantly,  the  models 
predict  moduli  which  have  the  same  frequency  dependence  as  is  observed 
in  frequency-dependent  moduli  in  typical  viscoelastic  materials.* 

The  outstanding  question  is  whether  or  not  the  parameters  of  the  RT 
and  RTG  models  can  be  chosen  to  describe  accurately  the  properties  of 
a  particular  viscoelastic  material.  The  construction  of  the  RT  model 
for  the  elastomer  3M-467  is  the  topic  of  the  following  section. 


*In  addition,  a  generalized  derivative  constitutive  relation  occurs 
in  Newtonian,  viscous  fluids  as  demonstrated  in  Appendix  B. 
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SECTION  V 

THE  RT  MODEL  FOR  THE  ELASTOMER  3M-467* 

The  first  material  examined  for  the  possible  application  of  generalized 
derivative  constitutive  relations  was  the  adhesive  tape  3M-467.  The  tape 
was  chosen  as  a  prime  candidate  because  of  its  viscoelastic  mechanical 
properties,  its  linear  response  in  shear  for  engineering  shear  strains 
up  to  1,  its  growing  applications  in  mechanical  damping,  and  the  fact 
that  sufficient  data  on  its  mechanical  properties  were  available. 

The  proposed  uniaxial  shear  RT  model  for  3M-467  is 

•  2  (-0  *  "/1>  Sm(t)  •  "  *  "  <"> 


where 

=  1.0  lb/in2 


(78) 


=  7.3  lb-sec' 56/in2 


(79) 


and 


(80) 


The  parameters  of  the  model,  p  ,  p-|  ,  and  a-j ,  are  chosen  so  the 
sinusoidal,  steady-state  response  of  the  model  closely  approximates  the 
sinusoidal,  steady-state  response  of  the  material  observed  experimentally. 
For  sinusoidal  strain 


c  (t) 
mnv 


e  iwt 
mnoe 

0. 


t  >  0. 
t  <  0 . 


(81) 


*3M-467  is  an  adhesive  produced  by  the  Minnesota  Mining  and  Manufacturing 
Co.,  Inc.,  Minneapolis,  Minnesota. 
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the  RT  model  generates  stresses  of  the  form 

'  2  ("o  *  “l11”5  £mn  elWt 

t  >  >0 


(82) 


for  t  large  enough  for  the  transients  to  have  died  out.  The  frequency- 
dependent  shear  modulus,  u(oj),  is  seen  to  be 

p(w)  =  u0  +  u^(iw)  (83) 

Figure  2  displays  the  good  agreement  between  the  experimentally  observed 
mechanical  properties  of  3M-467  at  75°F  and  the  RT  model  using  the  values 
of  the  parameters  given  above.* 

The  parameters  of  the  RT  model  were  determined  in  an  iterative 
manner.  Initial  guesses  of  the  parameters  were  made,  and  the  resulting 
frequency  dependent  shear  modulus  was  compared  to  the  observed  modulus. 
Successive  guesses  of  the  parameters  were  made  to  match  the  slopes  and 
asymptotes  of  the  model  to  those  of  the  observed  properties  until  an 
acceptable  fit  was  obtained. 

Although  the  parameters  of  the  RT  model  are  based  on  the  sinusoidal 
response  of  the  material  at  75°F,  the  model  can  be  used  for  non-periodic 
strain  histories.  To  demonstrate  the  ability  of  the  model  to  portray 
accurately  the  behavior  of  the  material  when  undergoing  non-periodic 
motion,  the  response  of  the  material  as  predicted  by  the  RT  model  is 
compared  to  the  experimentally  observed  response  of  the  material  at  75°F. 

In  particular,  the  behavior  of  3M-467  was  observed  when  the  material 
was  used  as  a  viscoelastic  spring  in  a  simple  oscillator  undergoing 


♦The  mechanical  properties  of  3M-467  were  provided  by  the  U.S.  Air  Force 
Materials  Laboratory,  Wright-Patterson  Air  Force  Base,  Ohio. 
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non-periodk  notion.*  The  viscoelastic  spring  was  two  pads  of  3M-467 
that  underwent  shear  strain  during  the  motion  of  the  oscillator.  Each 
pad  of  3M-467  was  made  by  laminating  2  mil  layers  of  3M-467.  Air 
entrapped  between  the  layers  of  the  pad  was  removed  by  pressing  the 
layers  together  with  a  5  lb.  weight  for  48  hours. 

The  other  two  components  of  the  oscillator  are  the  inass  and  the 
support  structure  (Figure  3).  The  mass  for  the  oscillator  is  a  metal 
cube  sandwiched  between  the  two  viscoelastic  pads  (Figure  4).  Each  pad 
is  attached  to  an  aluminum  brace.  Both  braces  are  glued  to  an  aluminum 
base  which,  in  turn,  is  glued  to  a  steel  foundation  as  shown  in  Figure  3. 

The  specific  objective  of  the  experiment  was  to  determine  the 
acceleration  transfer  function  of  the  oscillator.  The  acceleration 
transfer  function  is  the  ratio  of  the  transform  of  the  acceleration 
time  history  to  the  transform  of  the  input  force  history.  The  force 
time  history,  measured  by  a  Wilcoxon  Z-ll  impedance  head,  was  sampled 
at  2  x  104  measurements  per  second  and  all  frequencies  above  8  x  10^  Hz 
are  filtered  out.  The  mass  of  the  oscillator  was  tapped  with  a 
Wilcoxon  Z-ll  impedance  head  to  produce  impulsive  loading.  The  force 
time  history  measured  by  the  impedance  head  and  the  resulting  acceleration 
time  history,  measured  by  an  Endevco  accelerometer,  Model  2217,  were 

4 

also  sampled  at  2  x  10  measurements  per  second  where,  as  before,  all 

3 

frequencies  above  8  x  10  Hz  were  filtered  out.  The  transforms  of  the 
time  histories  were  calculated  using  the  "fast  Fourier  transform" 
routines  of  the  Hewlett  Packard  System  54 5 1 B . 

This  experimentally  determined  transfer  function  is  compared  to  the 
analytically  predicted  transfer  function  based  on  the  equations  of  motion 
of  the  oscillator  and  the  RT  model  for  the  viscoelastic  pads.  The  force- 
displacement  relation  for  the  two  pads  based  on  the  RT  model  (Equation  77) 
is 

fP(t)  =  x  *>0 +  x(o  w 

*A  schematic  of  the  viscoelastic  spring  appears  in  Figure  4. 
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Schematic  of  the  Oscillator  and  its  Support  Structure 
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where  A  is  an  area  of  contact  with  the  mass  for  each  pad,  6  is  the 
thickness  of  the  pad,  f  (t)  is  the  total  force  acting  on  the  faces  of 
the  pads  in  contact  with  the  mass,  and  x(t)  is  the  displacement  of  the 
face  of  the  pad  in  contact  with  the  mass.  This  force  displacement 
relation  is  based  on  the  assumption  that  the  displacement  in  the  pad 
varies  linearly  between  the  support  wall  and  the  mass  of  the  oscillator.* 


The  resulting  equation  of  motion  for  the  oscillator  is 

f(t)  =  mx(t)  *  (u0  +  x  ( t )  (85) 

Taking  the  Fourier  transform  of  the  equations  of  motion  and  determining 
the  acceleration  transfer  function  produces 


(ioo)  X(oi) 

“TO - 


2A 

6 


(u. 


(i<a) 


(iw) 


(86) 


A  comparison  of  the  experimentally  determined  and  analytically 
predicted  transfer  functions  for  five  oscillators  with  various  masses 
and  viscoelastic  spring  stiffnesses  is  presented  in  Figures  5  through  14.** 
Each  transfer  function  is  displayed  in  terms  of  its  magnitude  and  phase. 

The  agreement  between  the  observed  and  predicted  transfer  functions  is 
very  good. 

For  comparison,  the  calculated  transfer  function  based  on  a  Voigt 
viscoelastic  model  of  3M-467 

*  2  ("oc»n(t)  *  "l£mn(t))  •  "  *  n  (87) 


*A  finite  element  analysis  of  the  viscoelastic  pad  verifies  this 
assumption  to  be  valid  for  the  frequency  range  of  the  tests,  0  to 

5  x  103  Hz. 

**The  relevant  parameters  for  each  of  the  five  oscillators  are  given  in 
Table  1. 
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Includes  the  mass  of  the  accelerometer. 


(^as-qD/u.i  OL  -  9pnqLu6ew 
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f  the  Transfer  Function  for  Case  #1 


-he  Transfer  Function  for  Case  #1 


le  Magnitude  of  the  Transfer  Function  for  Case  #2 
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le  Magnitude  of  the  Transfer  Function  for  Case  #3 
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The  Phase  of  the  Transfer  Function  for  Case  # 3 


of  the  Transfer  Function  for  Case  H 


The  Phase  of  the  Transfer  Function  for  Case  *4 
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he  Magnitude  of  the  Transfer  Function  for  Case  #5 


Figure  14.  The  Phase  of  the  Transfer  Function  for  Case  #5 
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is  also  given  in  Figures  5  through  14.  The  parameters  of  the  Voigt  model, 

3 

.•o  and  |i,,  are  chosen  to  match  the  properties  of  3M-467  at  10  Hz. 


=  630  lb/in2 

(83) 

.113  lb-sec/in2 

(89) 

Note  that  for  some  of  the  oscillators,  the  Voigt  model  and  the  PT 

model  both  generate  transfer  functions  that  agree  reasonably  we1!  wii*> 

the  observed  transfer  functions.  However,  for  those  oscillators  having 

the  peak  magnitude  of  acceleration  response  at  higher  frequencies, 

(Figure  13  for  example)  the  transfer  function  based  on  the  RT  model  is 

clearly  in  better  agreement  with  the  measured  transfer  function  than  the 

transfer  function  based  on  the  Voigt  model.  In  addition,  the  phase  of  the 

observed  transfer  functions  is  consistently  modeled  more  accurately  by 

the  phase  of  the  transfer  functions  calculated  using  the  RT  model.  These 

results  follow  directly  from  the  fact  that  the  RT  model  accounts  for  the 

observed  properties  of  3M-467  over  the  entire  frequency  range  of  interest, 

2  3 

10  Hz  to  5  x  10  ,  whereas  the  Voigt  model  accounts  for  the  observed 

3 

properties  of  3M-467  only  in  the  neighborhood  of  10  Hz.  This  is 
clearly  seen  by  comparing  Figures  2  and  15. 

If  one  attempts  to  duplicate  the  results  presented  here,  one  should 
be  aware  that  the  mechanical  properties  of  3M-467  are  strongly  dependent 
on  the  water  present  in  the  material.  Figure  16  shows  the  variation  of 
the  real  part  of  the  modulus  with  relative  humidity.  Changes  in  the 
imaginary  part  of  the  modulus  with  relative  humidity  are  roughly 
proportional  to  changes  in  the  real  part.  Hence,  the  loss  factor,  the 
ratio  of  imaginary  part  to  real  part  of  the  modulus,  is  relatively 
insensitive  to  changes  in  relative  humidity,  as  seen  in  Figure  17. 

The  pads  of  3M-467  used  in  this  experiment  were  fabricated  under 
conditions  of  40%  relative  humidity  at  room  temperature.  However,  the 
pads  were  kept  covered  during  the  time  between  fabrication  and  installation 
into  the  test  setup. 
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Figure  17.  The  Variations  of  Loss  Factor  of  3M-467  with  Changes  in 
Relative  Humidity 
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1  „  summary ,  the  Rf  model  for  3M-467  is  capable  of  accurately 
predict iriy  the  non-per iodic  response  of  the  material  over  several  decades 
of  frequency,  and  is  superior  to  a  Voigt  model  of  the  material, 
consequently,  the  Ri  model  having  parameters  based  on  the  sinusoidal 
steady  motion  of  the  material  at  numerous  frequencies  is  capable  of 
predicting  the  response  of  the  material  to  impulse-like,  short  duration 
loading.  Therefore,  one  can  conclude  that  the  RT  model  can  accurately 
predict  the  general  response  of  the  material  within  the  frequency  range 
of  the  model. 
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SECTION  VI 

A  FINITE  ELEMENT  FORMULATION 
OF  THE  EQUATIONS  OF  MOTION 

Having  established  that  a  very  elementary  formulation  of  the  equation 
of  motion  for  an  elastomer-damped  oscillator  produced  excellent  agreement 
with  experimental  observation,  it  is  appropriate  at  this  point  to  put 
forward  the  tools  required  in  the  analysis  of  more  complicated  structures 
of  engineering  interest.  In  particular,  the  development  focuses  on  the 
analysis  of  structures  having  both  elastic  and  viscoelastic  components. 

A  continuum  formulation  of  the  equations  of  motion  for  such  structures 
is  impractical  because  of  the  resulting  complexity  of  the  formulation 
for  most  structures  with  complex  geometry  and  varying  material  properties. 
As  a  result,  a  finite  element  formulation  of  the  equations  of  motion  is 
adopted. 

The  cornerstone  of  the  finite  element  approach  is  the  construction 
of  the  stiffness  matrices  for  each  of  the  finite  elements  in  the  structure. 
The  stiffness  matrices  for  the  elastic  finite  elements  of  structure  are 
constructed  in  the  normal  fashion  using  assumed  displacement  methods  or 
assumed  stress  methods,  etc. 

The  formulation  of  the  stiffness  matrices  for  the  finite  elements  in 
the  viscoelastic  components,  however,  is  limited  to  those  methods  that 
do  not  constrain  the  stresses  in  each  finite  element  to  be  in  equilibrium 
with  the  forces  at  the  nodes  the  element.  The  assumed  stress  method,  in 
particular,  is  based  on  this  constraint  (Reference  24).  As  a  result, 
the  time  dependence  of  the  stresses  in  the  element  is  predicated  on  time 
dependence  of  the  nodal  forces.  However,  this  contradicts  the  fundamental 
nature  of  the  generalized  derivative  models  in  which  the  time  dependence 
of  the  stresses  is  predicated  on  the  time  dependence  of  the  strain 
histories. 
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c ,oc::t  i,  .  assumed  displacement  method  is  adopted  to 
,  . .  ;ii  I  a !  ia  tnc-  .t  .  t'tiits:  matrix  of  the  viscoelastic  finite  element 
(!'  ference  ^r>) .  m  an  assumed  displacement  element,  the  displacements 
v.iuiiri  the  element  are  assumed  functions  of  the  nodal  displacements. 

The  stilt i,t  mat i  ix  of  the  viscoelastic  finite  element  is  constructed 
I. .  I i 1 1 1  the  ei  -jitii.  viscoelastic  correspondence  principle.  The  stiffness 
matrix  is  first  formulated  as  though  the  material  were  elastic.  The 
.tiffness  matrix  is  then  separated  into  two  matrices,  one  matrix  containing 
those  elements  proportional  to  the  elastic  constant  A,  and  the  other 
matrix  containing  those  elements  proportional  to  the  elastic  constant  p. 

[K  ]  =  A[K‘]  +  p[K"]  (90) 

“  6 

At.  this  point  the  transforms  of  the  moduli,  u*(s)  and  A*(s),  from  either 
i: ht:  RT  or  the  RTG  models,  are  substituted  in  place  of  the  elastic  constants, 
and  >.  The  tesuli  is  the  viscostiffness  matrix  of  the  finite  element, 

I  i't  (  *  )  J 

IK,  Ul]  -  X*(s)[K']  +  u*(s)[K'']  (91) 

v  C  t: 

. ri. 1  •  i  f  i  .ill  a .  the  .  isv  /stiffness  matrix  for  a  finite  element  in  which  the 
hit)  const  i  hi  live  relation  is  used  to  model  the  material  is 


i  1  e  1  1  1 
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The  viscostiffness  matrix  of  the  finite  element  relates  the  nodal 
forces,  1 F(s) I ,  and  nodal  displacements,  { X ( s ) } ,  as  shown  below 

{F(s)  }  =  [Ke(s)]  {X  (s )  }  (93) 

and  the  viscostiffness  matrix  of  a  viscoelastic  structural  component  is 
constructed  from  the  viscostiffness  matrices  of  the  elements  within  the 
component  in  the  normal  manner.  The  viscostiffness  matrices  of  the  R 
viscoelastic  structural  components 

(F(s)}r  =  [K(s)]r  {X(s))r  r  -  1,2,3,...  ,R  (94) 

and  the  stiffness  matrices  of  the  Q  elastic  structural  components 

(F(s)}q  =  [  K]q  {X  (s )  }q  q  =  1,2,3,. ,.,Q  (95) 

are  used  to  construct  the  stiffness  matrix  of  the  total  structure,  again 
in  the  normal  manner. 

The  stiffness  matrix  of  the  total  structures,  [K(s)],  and  the  mass 
matrix  of  the  total  structure  are  now  used  to  construct  the  Laplace 
transform  of  the  equations  of  motion  of  the  structure 

s 2 [M]  {X (s ) }  +  [ K (s ) ]  {X (s ) }  =  (F  (s)  }  (96) 

Since  some  of  the  elements  in  [K(s)]  are  functions  of  s,  decoupling  the 
equations  of  motion,  (Equation  96)  to  obtain  solutions  is  more  complicated 
than  decoupling  the  equations  of  motion  of  a  completely  elastic  structure 
where  the  stiffness  matrix  has  constant  elements.  Finding  solutions  to 
Equation  96  is  the  topic  of  the  next  section. 
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SECTION  VII 

THE  SOLUTION  OF  THE  DISCRETE  EQUATIONS  OF  MOTION 

The  task  at  hand  is  the  solution  of  the  equations  of  motion  which 
resulted  from  the  finite  element  formulation  (Equation  96).  A  form  of 
modal  analysis  is  adopted  where  the  mode  shapes,  eigenvectors,  of  the 
equations  of  motion  are  used  to  construct  an  orthogonal  transformation 
of  the  variables  that  decouple  the  equations  of  motion.  The  decoupled 
equations  of  motion  are  then  used  to  determine  the  components  of  the 
structure's  response  and  a  general  form  of  the  solution  to  the  equations 
of  motion  is  derived. 

Throughout  this  development,  the  viscoelastic  components  of  the 
structure  are  described  by  their  respective  RTG  models.  Since  the  RT 
model  is  a  simplified  version  of  the  RTG  model,  the  method  of  solving 
the  equations  of  motion  of  a  structure  with  viscoelastic  components 
described  by  RT  models  will  be  seen  to  be  a  special  case  of  the  following 
solution  technique. 

The  reason  for  developing  a  special  solution  technique  for  the 
equations  of  motion 

s 2  [M]  { X (s )  }  +  [K  (s )  ]  {X  (s)  }  =  { F  C s ) }  (97) 

is  that  the  normal  method  of  decoupling  the  equations  of  motion,  using 
modal  analysis  to  construct  an  orthogonal  transformation  that  diagonalizes 
the  mass  and  stiffness  matrices,  is  not  applicable  because  the  stiffness 
matrix  of  the  structure,  [K(s)],  contains  terms  that  are  dependent  on 
the  Laplace  parameter,  s. 

The  method  of  solution  for  Equation  97  is  an  extension  of  the  method 
proposed  by  Foss  to  decouple  the  equations  of  motion  for  a  structure  with  . 

non-proportional  viscous  damping  (Reference  26).  j 

i 

1 

Cm]  +  [C]  (x  ct ) }  +  [K]  { x  (t )  }  =  { £  (t ) }  (98) 
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Non-proportional  damping  occurs  when  the  damping  matrix  [C]  is  not  a 
linear  combination  of  the  mass  and  stiffness  matrices  of  the  structure. 

[C]  /  ajD-1]  +  a2[K]  (99) 


At  present  there  is  no  general  method  of  constructing  an  orthogonal 
transformation  for  three,  real,  square,  symmetric  matrices  when  each 
of  the  matrices  is  not  a  linear  combination  of  the  remaining  two. 
Consequently,  Foss  posed  the  equations  of  motion  for  non-proportional 
damping  in  terms  of  two  real,  symmetric  matrices. 


d 

ar 


'  O  I  M  ' 

X 

+ 

■-M  |  O  ‘ 

.  M  1  C  . 

.  X 

-OIK. 

X 

X 


0 


f  (t) 


(100) 


1 


The  lower  set  of  the  partitioned  matrix  equations  is  the  equations  of 
motion  of  the  structure  and  the  upper  set  of  matrix  equations  is  satisfied 
identically.  The  equations  of  motion  as  posed  in  Equation  100  are  readily 
decoupled  and  solved. 

To  solve  the  equations  of  motion  of  the  structure  containing  elastic 
and  viscoelastic  components  (Equation  96)  the  equations  are  posed  in  terms 
of  two  real,  square,  symmetric  matrices.  To  begin,  one  multiplies  the 
equations  of  motion  by  each  distinct  term  appearing  in  the  denominators 
of  the  elements  of  the  stiffness  matrix,  [K(s)j, 

[D(s)s?[M]  +  D(s)[K(s)]]  {X  (s ) }  =  D(s)  (F(s)}  (101) 

where 

N  K°  6nV  N  ?n  L 

DCs)  =  n  (1  *  l  a  ,s  nk)  •  n  (1  +  l  b  s  np)  (102) 
n=l  k=l  nK  n=l  p=l  np 
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assuming  that  there  are  N  different  viscoelastic  materials  in  the 
structure.  Multiplying  D(s)  times  [K(s)L  the  stiffness  matrix,  produces 
a  matrix,  [Kp(s)],  that  has  no  terms  in  s  appearing  in  the  denominators 
of  its  elements. 

D(s)  [K(s)]  =  [Kd(s)]  (103) 

In  fact,  all  of  the  elements  of  [K^(s)]  are  constant  terms  plus  terms 
containing  s  raised  to  real,  positive  powers.  Also  note  that  the  matrix 
D(s)s  [M]  has  elements  which  are  sums  of  terms  containing  s  raised  to 
real,  positive  powers. 

The  equations  of  motion  are  now  expressed  as 

[  Z  (s )  ]  {X  (s ) )  =  D(s){F(s)}  (104) 

where 

[  Z  (s )  ]  =  [D(s)s2[M]  -  [  Kj)  (s )  ]  ]  (105) 

At  this  point  in  the  development,  the  real,  positive  exponents  of  s 
appearing  in  the  matrix  [Z(s)]  are  taken  to  be  rational  as  well.  Had 
any  of  the  exponents  been  initially  irrational,  they  are  replaced  by  their 
rational  approximations  to  as  many  significant  digits  as  desired.  Since 
all  of  the  exponents  in  [Z ( s ) ]  are  rational,  the  matrix  may  be  expressed  as 

[Z(S)]  =  [  [Mj  i  c.sJ/m  *  l  [K  JsVm  ]  (106) 

j  =2m^  *=1  1 

where  m  is  the  smallest  common  denominator  of  the  exponents  of  s  in 
[Z(s)J  and 

J  3/ 

s 2D (s )  =  l  c.s  /m  (107) 

j  =  2m  ^ 
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[Kd(s)] 


L  Vm 
l  C  k  3  s  m 
4-0  l 


(108) 


where  [K?]  is  symmetric  and  some  of  the  [K^]  and  appearing  above 
may  be  zero. 


Using  Equation  106,  the  equations  of  motion  of  the  structure  become 


[  [M]  j  c.sj/m  *  l  [K.]si/m]  {X  (s ) }  =  D(s){F(s)}  (109) 
j  =  2m  ■*  t  =  0 


and  expressed  in  terms  of  one  index  of  summation,  they  are 


J  j  / 

l  [  CM]ci  +  [  K  -  ]  ]  s  m  { X  (s )  }  =  D(s){F(s)J  (110) 

i  =  0  J  J 


or 


J  j  / 

l  [A  -  ]  s  m  {X  (s ) }  -  D(s)(F(s)}  (111) 

j  =  0  J 


where 


[Aj  ]  =  [  [M]cj  +  [K.J  ]  (112) 

and  again  recognizing  that  some  of  the  c.  and  [K.]  are  zero. 

J  J 

The  equations  of  motion  as  given  in  Equation  111  are  now  posed 
in  terms  of  two  real,  square,  symmetric  matrices. 


* 

t 

-  ■ .  -1 

A 
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f  'n 

(0.) 

(0.) 

{0.} 

'V/  . 

=  *  .  -v  (113) 

(0.) 

{0.} 

D(s){F(s) } 

v  J 

Note  that  each  matrix  [A.],  j  =  0,1,2,...,J,  is  real,  square,  and 
symmetric  because  they  are  linear  combinations  of  [M]  and  [K.].  It 

J 

follows  that  the  two  matrices  containing  [A,]  in  Equation  113  are  also 

J 

real,  square,  and  symmetric.  Also  notice  that  the  lowest  set  of 
partitioned  matrix  equations  in  Equation  113  are  the  equations  of  motion 
of  the  structure  as  given  by  Equation  111,  and  that  all  of  the  upper  sets 
of  matrix  equations  in  Equation  113  are  satisfied  identically. 

The  equations  of  motion  as  posed  in  Equation  113,  referred  to  as  the 
expanded  equations  of  motion,  can  be  decoupled  using  an  orthogonal 
transformation.  The  general  form  of  the  expanded  equations  of  motion  is 

s?-[M]  {X  (s) }  +  [K]{X(s)J  =  (F(s ) }  (114) 

which  is  Equation  113  expressed  in  more  compact  notation.  The  orthogonal 
transformation  is  constructed  from  the  eigenvectors  associated  with  the 
eigenvalue  problem  for  the  expanded  equations  of  motion. 

Xn[M]{^n>  +  mUn>  =  {°->  015) 
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The  eigenvectors  [ij>nl  are  used  to  construct  the  orthogonal  transformation 
matrix  [<t>]  in  the  normal  fashion,  and  the  resulting  transformation  of 
variables  is 

(X  (s) }  =  [*]{a(s)}  (116) 

Substituting  this  transformation  into  Equation  14  and  premultiplying  the 
equation  by  [4>]^  produces 

s1/m[i]T[M][>]{a(s)}  +  [^[kHi]  a  (s )  =  I>]T{F  (s) }  (117) 

To  demonstrate  that  Equation  117  is,  in  fact,  the  decoupled  expanded 
equations  of  motion,  one  uses  the  fact  that  eigenvectors  of  the  expanded 
equations  of  motion  are  orthogonal  with  respect  to  [M]  and  [K]. 

{*j  )T[M3Un}  =  0.  j  t  n  (118) 

U  }T[K]Un}  =  0.  j  ^  n  (119) 

Equation  117  then  reduces  to 

i 

s'^N^afs)}  +  [^kr]{a(s)}  =  [$]T{F(s))  (120) 

where  ]  and  [ —  kn —  1  are  diagonal  matrices  of  the  modal  constants 

and  k^,  respectively. 

mn  =  UnlT[M]('n}  (121) 


Premultiplying  Equation  120  by  C  ' mn — ] “ ^  or  equivalently  ["•- — .]  yields 

|( 

s1/m[-i^]{a(s)}  +  [-_!!.}[  a  (s)>  =  ["^-0  [  If  (F  (s  )  }  ( 1 23) 

n  n 

i,L  i.L 

The  ratio  of  the  ntn  modal  parameters,  kn/mn,  is  minus  the  n1" 
eigenvalue  of  the  expanded  equations  of  motion. 
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Premultiplying  Equation  115  by  { 4>n }  produces 

Xn{  }T[M]  { 4>n }  +  Un}T[K]Un)  -  0 


X  m  +  k 
n  n  n 


from  which  Equation  124  follows.  As  a  result,  the  equations  of  motion 
can  be  further  simplified  to 

s  {a  (s )  }  -  [-A^]{a(s)J  =  [-'SpJ[4]T{F(s)}(127) 

n 

From  Equation  127  it  follows  that  the  expression  for  the  Laplace 
transform  of  the  n^  modal  coefficient,  an(s),  is 


an(s) 


{ 4>n  )T{F(s)  } 


This  expression  for  the  modal  coefficient  can  be  further  simplified  by 
noting  the  general  form  of  the  eigenvector  { 4>n }  associated  with  the 
expanded  equations  of  motion 

r  > 

xJ'2{<t>  } 

n  Tn 


n  n 


w 
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where  {>)>  f  and  A  are  the  solutions  of  the  eigenvalue  problem  associated 
with  the  original  equations  of  motion  in  the  form  of  Equation  111. 


{0.  }  n  =  1,2,3,  .  .  .  ,N  (130) 


The  general  form  of  the  ntn  eigenvector,  {6nl,  can  be  verified  by  direct 
substitution  in  Equation  115  when  [M]  and  [K]  are  expressed  in  terms  of 
the  matrices  [A^]  as  indicated  in  Equation  113.  The  lowest  set  of 
resulting  partitioned  matrix  equations  produces  Equation  130,  and  the 
upper  sets  of  the  partitioned  matrix  equations  are  satisfied  identically. 
Consequently,  the  numerator  of  the  n^  modal  coefficient,  (4>  l^f F(s) } ,  is 


Un}T{F(s)} 


r 

-j-1, 

T 

r 

xn 

{0.  } 

-J-2, 

Xn  Un} 

(0.  } 

X,T{V 

(0.  } 

%  *  %  # 

'V  %  'X*  ,  % 


« 

• 

A  { $  } 

n  n 

{0.  } 

V*n} 

{0.,} 

(*n} 

D(s)  {F  (s ) } 

L  n  J 

^  J 

(131) 


or 


{^n}T{F(s)}  =  Un}T{F(s)}DCs) 


(132) 
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and  the  n^  modal  coefficient  reduces  to 


an(s) 


Un)‘{F(s))D(s) 

^  /m  ~  " 
mri  s  ■An) 


(133) 


The  Laplace  transform  of  the  displacement  response  of  the  structures 
follows  from  Equation  116  and  takes  the  form 


N 

(X  (s )  }  *  V  an(s){*n)  (134) 


or 


N  U  }  (F(s)}D(s) 
(X  (s )  }  l  — - —  1-— - 


n=l 


m  (s  '  m- X  ) 
n  n' 


(135) 


where  N  is  the  order  of  the  matrices  [M]  and  [K]  in  the  expanded  equations 
of  motion. 

The  order  of  the  expanded  equations,  N,  can  be  very  large.  From 
Equation  113  it  is  clear  that  the  order  of  the  expanded  equations  of 
motion  is  equal  to  6,  the  order  of  the  matrices  [Aj],  times  J  where  J  is 
defined  in  Equation  107. 

N  =  6  •  J  (136) 

From  Equation  107,  it  is  clear  that  J/  is  the  largest  exponent  in  the 

r\  III 

expression  s^D(s)  which  is  2  +  8,  where  B  is  the  largest  exponent  in 
D(s).  Therefore, 

J  «  m  (2  +  e)  (137) 

and  the  order  of  the  expanded  equations  is  seen  to  be 

N  =  6m  (2+  $)  (138) 

Note  that  if  m,  the  smallest  common  denominator  of  the  rational  exponents 
of  s,  in  the  original  equations  of  motion  is  large,  the  order  of  the 
expanded  equations  is  quite  large  for  a  structure  with  anything  more  than 
a  very  modest  number  of  degrees  of  freedom.  However,  the  solutions  to  the 
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expanded  equations  of  motion  can  be  obtained  by  using  numerical  methods 
that  do  not  involve  the  manipulation  of  the  expanded  equations,  as  will 
be  shown  in  the  following  section. 

Before  proceeding,  it  should  be  pointed  out  that  the  equations  of 
motion  of  a  structure  containing  both  elastic  and  viscoelastic  components 
can  be  solved  given  that  a  finite  element  formulation  of  the  equations 
of  motion  is  possible  and  that  an  RT  or  RT6  model  exists  for  each 
viscoelastic  material  in  the  structure.  The  general  form  of  the  solution 
technique  for  the  equations  of  motion  containing  only  RT  models  for  the 
viscoelastic  components  is  identical  to  the  above  development  except  that 
D(s)  is  set  equal  to  one. 
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SECTION  VIII 

CALCULATING  THE  LAPLACE  TRANSFORM 
OF  THE  STRUCTURAL  RESPONSE 


At  this  point,  it  is  clear  that  solutions  to  the  equations  of  motion 
for  the  structure  containing  elastic  and  viscoelastic  components 


(X(s) } 


„  Un>  (F(s)}D(s) 

1  r  T  - 

'  m  -  \  1 
m^Cs  nJ 


N 

l 

n=l 


(135) 


are  difficult  to  calculate  from  the  expanded  equations  of  motion,  because 
of  the  large  order  of  the  matrices  in  the  expanded  equations.  As  a 
result,  an  alternative  method  of  obtaining  the  solutions  is  required  if 
the  finite  element  formulation  of  the  equations  of  motion  using 
generalized  derivative  models  is  to  be  a  useful  tool  to  the  engineer. 


The  alternative  method  adopted  here  is  a  combination  of  iterative 
schemes  used  to  obtain  the  eigenvalues,  An>  and  eigenvectors,  {<t>n), 
associated  with  the  original  equations  of  motion. 


J  -i 

[  l  [A.]aJ]U}  =  {0.}  ,  n  =  1,2,3,.  ..,N  (130) 

j  =  0  J  n  n 


Recall  that  the  number  of  distinct  homogeneous  solutions  to  the  equations 
of  motion,  N,  is  dependent  on  the  smallest  common  denominator  of  the 
exponents  in  the  equations  of  motion,  m;  the  largest  exponent  in  the 
product  of  the  denominators  terms  of  the  Laplace  transform  of  the  RTG 
models  in  the  equations  of  motion,  6;  and  the  number  of  degrees  of 
freedom  of  the  structure,  <5. 

N  =  6m(2  +  f?)  =  2m6  +  6m6  (138) 


63 


AFML-TR-79-41 03 


Using  an  iterative  scheme  based  on  the  homogeneous  form  of  the 
equations  of  motion  in  Equation  96,  26m  of  the  homogeneous  solutions 
are  obtained. 

Upm[M]  +  [K(Xn)]Un}]  =  {0.}  (139) 


The  iteration  for  the  solutions  centers  on  calculating  successive 
estimates  of  >n  and  {<f.  j  using  the  scheme 

[-2m(P+l)[M]  +  CK(^P))]]{«n)(P+1)  =  (0.)  (140) 

A^  is  the  pth  estimate  of  A  .  >(2m(P+1)  is  the  (p  +  l)*"*1  estimate  of 

A^m  and  is  the  (p  +  l)*'*1  estimate  of  { <pn } . 


Given  a  value  of  A^,  the  (p  +  l)th  estimates  of  A^m  and  { <j3n >  may  be 
calculated  using  matrix  iteration  or  any  other  method  that  is  appropriate 
to  obtain  the  solution  to  the  eigenvalue  problem 

=  uU}  (141) 


where  B  is  complex.  Note  that  Equation  140  can  be  expressed  as 


[K(x(P))r1[M]Un}(Ptl) 


1  ,(P+1) 

piTFHT  {V 

n 


(142) 


which  is  of  the  same  general  form  as  Equation  141. 


The  iterative  scheme  as  it  appears  in  Equation  142  is  only  useful  in 
obtaining  the  eigenvalue  with  the  smallest  magnitude,  A-j ,  and  the 
associated  eigenvector,  To  obtain  the  other  eigenvalues  and 

eigenvectors,  the  iteration  scheme  is  modified  to  allow  the  scheme 
to  converge  on  the  larger  eigenvalues  and  associated  eigenvectors. 

For  instance,  the  scheme  used  to  obtain  the  eigenvalue,  A.  ,  and 

±.  u  L 

the  Ltn  eigenvector,  (^},  is  (Reference  27) 
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[K(^P))]_1 


L 

I 

1=1 


H, 


« 2mTFj 


[m](.l)(p+1) 


1 


~2m(P+l ) 
AL 


<*L} 


(P+D 


(143) 


The  terms  in  the  summation  on  the  index  i  subtract  the  components  of 
the  first  L  -  1  eigenvectors,  associated  with  eigenvalue  problem, 

[K(x[p))]'1[M]{^}(p)  =  -  -liTTPy  {V(P)*=l,2,3,...>L-l(144) 

At 

from  the  successive  approximations  of  {$L}^P+^  calculated  when  matrix 
iteration  is  used  to  solve  Equation  143. 


Given  that  matrix  iteration  has  successfully  produced  A.  and 

( P+l  )  L 

{^}V  ,  one  can  take  advantage  of  the  orthogonality  of  the  mode  shapes 

T 

Wpl(P)  [M]{*.  >(P+1)  =  0.  1-1,2 . L-l  (145) 


to  demonstrate  that 


(P) 


1  ^  '  <»•> 


(146) 


and  Equation  143  reduces  to 

[K(»(P,}]'1[M]{*l}(P*1)  =  -  -2m(P,l)  <«Ll(P*1)  (’47) 


Thus,  {<bL } ^ P+1  ^  and  X^171^ P+1 )  are  in  fact  the  (p  +  -jjth  approximat:ion  of 
the  solution  to  the  equations  of  motion. 
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On  the  other  hand,  and  ^  =  1.2,. ...L-l,  are  not  in 

any  sense  approximations  of  the  solutions  to  the  equations  of  motion. 
However,  they  are  the  first  L-l  eigenvalues  and  eigenvectors  associated 

with  Equation  147.  This  can  be  seen  by  comparing  Equation  144  with 

(Pi  (Pi 

Equation  147.  {i^)v  '  and  '  can  be  calculated  using  matrix  iteration. 


Notice  that  in  continuing  the  iterative  processes  in  Equation  140  or 

143  the  (P  +  2)tfl  approximations  of  the  solution,  x2m(P+2)  and  {4  }^P+^, 

are  based  on  the  value  of  Xn'  .  However,  the  function  Z  m  has  2m 

branches.  Given  a  value  of  one  can  calculate  2m  values  of 

A  one  value  for  each  branch  of  Z^2m 

n 

So,  when  using  Equations  140  and  143  to  obtain  solutions  of  the 

1  /  pm 

equations  of  motion,  it  is  necessary  to  choose  one  branch  of  Z  to 

calculate  the  6  eigenvalues  and  5  eigenvectors.  Then  another  branch  of 
1/9 

Z  is  chosen  and  <5  other  solutions  are  obtained.  This  process  is 
continued  until  all  2m  branches  have  each  been  used  to  calculate  6 
solutions  producing  a  total  of  2m6  homogeneous  solutions  to  the  equations 
of  motion.  The  general  form  of  the  equations  that  has  these  2m6 
homogeneous  solutions  is 


+  (k)} 


j  =  1 , 2 ,  .  .  ,  ,  6 
{0.}  (148) 

k=l  ,  2 , . . . , 2m 


where  the  subscript  k  denotes  the  branch  of  Z 


l/2m 


on  which  the  relation 


is  valid. 


,~2m  ^  2ni 

1  J  00J 


Xj  (k) 


(149) 
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The  remaining  0mS  of  the  N  homogeneous  solutions  to  the  equations 
of  motion  are  determined  using  Equation  138,  an  iterative  scheme  based 
on  the  homogeneous  form  of  Equation  109.* 


[[M]  l 


j  =  2m 


c.xJ 
J  n 


I  [K £]x£]{*  } 

Jt=0  £  K  n 


{0.} 


(150) 


The  iterative  scheme  is 


J-i-  - i CP) 

+  y  xc.xJ  ) 

•  ->  J  n 
j  =  2m  J 


*  I  CK  ]X^(P)]{4  }(P+1)  =  (0.  }  (151) 

o=n  »  n  n 


or,  in  more  compact  form 
[[M]Q(Xn(P+1)Xn(P)) 


[KD(X^)]]Un>^^ 


(0.  }  (152) 


where,  as  before,  X^  is  the  Pth  estimate  of  the  eigenvalue,  X^P+1^ 
is  the  (P  +  l)th  estimate  of  the  eigenvalue  and  {<j>n)^P+^  is  the  (P+l)tfl 
estimate  of  the  eigenvector.  Given  a  value  for  X^p',  one  can  calculate 
1$  }^P+^  and  Q(^P+^,  Xn^)  using  matrix  iteration  or  any  other  method 
suitable  for  the  solution  of  Equation  141.  Equation  152,  expressed  in 
the  form  of  Equation  141,  is 


1 


[M]  { <f>  } 


(P+1) 


n  ■ 


{ ) 


(P+1) 


(153) 


♦Note  that,  if  the  equations  of  motion  contain  only  RT  models  for  the 
viscoelastic  materials,  that  D(s)  is  one  and  6,  the  largest  exponent 
of  s  in  D(s),  is  zero.  Hence,  the  total  number  of  homogeneous  solution, 
N,  is  2rri'5  as  seen  by  Equation  138.  Since  Equations  140  and  143 
provided  2m<S  solutions,  one  can  return  to  Equation  136  and  calculate 
structural  responses. 
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At  this  point,  the  numerical  value  of  Q(A 


r(P+l )  ,(P) 


is  used  to  calculate  A 


(P+1) 


from  Equation  153 


with  the  relation 


Qc;ri)-‘T)>  - 


J  n  n 


(?) 


J 

l 

j  =  2m 


c.A^) 
3  n 


(154) 


which  follows  from  Equation  151  and  152.  However,  this  method  of 
calculating  A^P+1 ^  assumes  that  >^p+^  are  both  on  the  principal  branch 

of  (Apm)  The  method  of  calculating  A^P+^  assuming  that  A^P+^  and 

A^P'  are  both  on  the  kth  branch  of  (A^m)^^m  ^ 


QCx 


: (P+i) 


n(k) 


A(P)  )  - 
Xn(k)J 


c  A J  * 

3  n  (k) 


J 

l 

3=1 

r^+ry 


i/ 


C  A~CP) 

.1  n(k) 


Bm 


6m 


,'(P+1) 

Xn(k) 


(155) 


where  the 


th 


branch  of  Z 


1/  Bm 


is  used  to 


calculate  A^j^ , 


The  resulting  form  of  the  iteration  process  is 


[KDUn(k)]]  1[M^{<J,n(k) 


(P+1) 


'  (PVT)  (Fr 
^un(k)/n(k) 


) 


Pn(k) 


(P+1) 


0  56) 


~  ~ (P+1 )  ~ ( P ) 

where  successive. estimates  of  A^kj  are  calculated  using  Q(A^kj  >An(k)^ 

from  Equation  156  and  then  using  Equation  155  to  calculate  A^j^.  When 

using  matrix  iteration  to  solve  for  Q(A^Pkj ^ ’^n(k)^  and  ^n(k)^  P+1  in 
Equation  156,  one  usually  obtains  only  the  Q  with  the  smallest  magnitude 
and  its  associated  eigenvector, 


To  obtain  successive  estimates  of  the  other  Q's  with  larger 
magnitudes  and  their  associated  eigenvectors,  the  iteration  scheme 
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•v/p+lj  ~(P) 

is  modified  as  before.  The  iterative  scheme  that  yields  Q(\(k)  ’AL(k) ^ 
and  is 


rx  fx(p)  n'1  +  y  _ rM]u  ifp+1) 

LVXn(k))J  L  *f.(W)  .OTT  LMJUL(k)} 

*-1  yCAt(k)  ’A£(k)J 


(P+1) 


^UL(k)  ,XL(k)J 


where  Q(A^  P*  j  ^,A^  P,L )  and  are  the  solutions  to 

tKD(;L(k))TltM](*t(l)i(p)  •-0,;TP«rrTO-Tl<’»oo)('>)  ow) 

Q(At(k) ,At(k) j 


(P+-]  \  /p) 

having  the  L  -  1  smallest  values  of  Q.  The  values  of  Q(A„  / ,  v /  ,A: , .  > 


and  (^£(1^)} yr}  Cdn  be  obtained  using  matrix  iteration. 


%  IKI  '  XI  XI 


Using  Equations  155,  156,  and  157,  66m  homogeneous  solutions  to  the 
equations  of  motion  are  obtained  for  each  branch  of  Z  Since  there  are 

gm  branches  of  Z^6"1,  the  iteration  process  should  produce  6m6  homogeneous 
solutions.  These  homogeneous  solutions  satisfy  the  homogeneous  equations 
of  motion. 

J  L  .  „  p-1,2,... ,6 


C[M]j2mCjXP(k)  +  JoCK*3W{W  =  {°'} 


p=l , 2 , . . . , 6  ^ 
k=l , 2 , . . . , 6m 


The  two  iterative  schemes  used  to  obtain  the  N  solutions  to  the 

homogeneous  equations  are  considered  to  have  converged  when  successive 

~/pi  ~(P+1) 

approximations  of  *n(k) (*„(£)  and  k)  )  are  approximately  the  same 
complex  number.  -(P+l) 
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The  general  convergence  criteria  of  the  iterative  schemes  are  considered 
beyond  the  scope  of  this  investigation;  however,  the  schemes  have  converged 
for  numerous  structures  considered  by  the  author.  Homogeneous  solutions 
for  an  example  problem,  calculated  using  the  iterative  scheme  given  above 
appear  in  Table  2. 


Note  that  all  of  the  parameters  of  the  Laplace  transform  of  the 
structure's  response  are  determined,  except  the  modal  constant  mn 
defined  by  Equation  121. 

mn  =  {^n)TtM]{^n}  (121) 


The  eigenvector  of  the  expanded  equations  of  motion,  },  can  be 
constructed  from  the  nLn  eigenvalue  and  associated  eigenvector  of  the 
original  equations  of  motion,  A  and  { 4>n } »  using  the  relation  given  in 
Equation  129.  This,  coupled  with  the  general  form  of  [M]  given  in 
Equation  113,  produces  an  expression  for  mn  which  takes  the  form 


f  <*. 


iTr  v  i 

i  '  L  J  n  *-  *-»  J  J  '  r  t 

j  =  i  " 


60) 


The  modal  constant  of  the  expanded  equations  of  motion,  mn,  can  be 
calculated  without  manipulating  the  expanded  equations  of  motion. 


In  conclusion,  all  of  the  parameters  in  the  general  form  of  the 
Laplace  transform  of  the  structure's  response  can  be  calculated  without 
manipulating  the  expanded  equations  of  motion,  given  that  the  iterative 
schemes  outlined  above  converge. 
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TABLE  2 

.  HOMOGENEOUS  SOLUTIONS  OBTAINED  USING 
THE  PROPOSED  ITERATIVE  SCHEMES  FOR  AN  EXAMPLE  PROBLEM 


[M] 


.3300  .0825 

.0825  .3300 


4.0  +  0,2  +  °-2s^  -2.0 

[K(s)]  = 

1.0  +  0.1sH 

-  2.0  4.0  +  0,01 

+  0 . 01 s  ^ 

1.0 

+  0 .  Is ^ 

Using  Equations  140 

and  143 

'1.000  +  i 

Al(l)  =  1-071768  +  i  1.099794 

1.040  +  i 

rl.000  +  i 

Xl(2)  =  1-  073498  +  1  1-028611  =1 

,1.001  +  i 

fl.000  +  i 

X1  4)  =  l-073498  '  1  1-028611  {<f>lf3)}  = 

ll.OOl  -  i 

fl.000  +  i 

Xl(5)  =  1.071768  -  i  1.099794  =  < 

L 1 . 040  -  i 

fl.000  +  i 

A2(1j  =  1.  581998  +  i  1.601989  t  <P  2  ( 1 )  }  =' 

,-.9713+  i 

I 

1.000  +  i 

*2C2)  =  - 1 .  5846  31  +  i  1.  547150  U2(2)}  =] 

.-1.002  +  i 

0.0 

0.0165 

0.0 

0.0238 

0.0 

0.0238 

0.0 

0.0165 

0.0 

0.0120 

0.0 

0.0240 
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SECTION  IX 

THE  EXISTENCE  OF  THE  STRUCTURAL 
RESPONSE  TO  IMPULSIVE  LOADING 


Of  particular  interest  at  this  juncture  is  whether  or  not  the  inverse 
transform  of  the  Laplace  transform  of  the  structure's  displacement 
response  for  impulsive  loading  exists.  In  fact,  the  inverse  transform 
always  exists  and  is  real,  continuous,  and  causal. 


To  demonstrate  this,  one  starts  with  the  general  form  of  the  transform 
of  the  structural  response. 


(X (s ) } 


¥  t  }T{ F  (s  )  }D  (s  ) 
n=l  17— -  {*n} 


m  (s 
n 


m 


-V 


(135) 


For  simultaneous  unit,  impulsive  loading  at  the  structure's  degrees  of 
freedom,  the  column  vector  of  applied  forces  is 

{ f  (t  )  }  =  6  (t )  { 1 .  }  (161) 

where  {1.}  is  a  column  vector  of  ones.  The  Laplace  transform  of  the 
column  vector  of  forces  is 

{ F  (s ) }  =  (1 .  }  (162) 

and  the  transform  of  the  response  for  the  impulsive  loading  is 


(X  (s ) } 


N  T 

"  U  )  (1.  }D(s) 

-—17--.-  -  <♦„> 


I 

n=l 


(163) 


The  inverse  transform  of  this  expression  always  exists,  which  follows 
from  a  theorem  on  the  existence  of  the  inverse  transform  (Reference  28). 
Paraphrasing  the  theorem  in  terms  of  the  notation  used  above;  it  states 
that  the  inverse  transform  of  { X( s ) >  exists  and  is  real,  continuous, 
and  causal  when 

1.  { X ( s ) }  is  analytic  for  Re[s]>0, 
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2.  : X ( s ) }  is  real  for  s  real  and  positive,  and 

3.  iX(s)i  is  order  s  Y,  where  y>l ,  for  s  large  in  the  right 
half  s  plane. 

;  X ( s ) }  is  analytic  for  Re[s]>0  when  the  branch  cut  of  s^™  is  chosen 
to  lie  along  the  negative,  real  axis  in  the  s  plane  and  the  poles  of 
{ X ( s ) },  which  occur  at 

s  =  Xm  (164) 

n  v 

and  do  not  appear  in  the  right  half  s  plane.* 

{ X ( s ) }  is  real  for  s  real  and  positive.  The  only  quantities  appearing 

in  Equation  163  that  can  be  complex  are  {<)>},  A  and  m  ,  because  D(s)  and 
^  n  n  n 

1  /  m 

s  are  real  for  s  real  and  positive.  When  {<Pn> ,  \  and  mn  are  complex, 
they  occur  in  conjugate  pairs.  Note  if  An  and  ($n)  are  a  homogeneous 
solution  to  the  equations  of  motion 

[  l  [A.]x-hu  }  =  {0.}  (130) 

j=l  3  n  n 

where  X  and  {*  }  are  complex,  that  their  conjugates  are  homogeneous 
solutions  to  the  equations  of  motion. 

[  I  [A.]Ij]{*  1  =  {0.}  (165) 

j=l  J  n  n 

Equation  165  follows  directly  from  the  complex  conjugate  of  Equation  130. 
Since  homogeneous  solutions  occur  in  conjugate  pairs,  the  modal  constant 
mn  occurs  in  conjugate  pairs. 

mn  =  {*n}Tti1j^n‘1CAj]]{^n}  (160) 


*A  pole  in  the  right  half  s  plane  indicates  that  an  RT  or  RTG  model  in 
the  equations  of  motion  characterizes  the  viscoelastic  material  as 
generating  energy  instead  of  dissipating  energy. 
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mn  '  {*n}T[  (166) 

It  follows  directly  that  when  the  terms  in  Equation  163  are  complex,  they 
occur  in  complex  conjugate  pairs  and  { X ( s ) }  is  real  for  s  real  and 
positive. 

{X(s)f  also  satisfies  the  third  and  last  condition  placed  on  the 

~  -2 

transform.  To  show  that  { X ( s ) }  is  order  s  for  s  large  in  the  right 
half  s  plane,  one  uses  the  transformed  equations  of  motion  as  they 
appear  in  Equation  97. 

s2[M]  {X  (s)  }  +  [K(s)]{X(s)l  =  {F  (s ) }  (97) 

The  transformed  equations  of  motion  for  simultaneous,  unit  impulsive 
loading  is 

[s 2  [M]  ♦  [K (s )  ] ]  (X  (s )  }  -  (1.1  (167) 

The  only  terms  in  [K(s)],  other  than  the  constant  terms,  are  those  terms 
proportional  to  p*(s)  and  A*(s)  from  the  RT  and/or  RTG  models  of  the 
viscoelastic  materials.  The  general  forms  of  p*(s)  and  X*(s)  appear  in 
Equation  62  for  the  RTG  model  and  Equation  63  for  the  RT  model .  As  a 
direct  result  of  the  general  forms  of  n*(s)  and  A*(s),  Equation  167 
reduces  to 

s 2  [M]  (X  (s )  }  =  (1.)  (168) 

for  s  large.  Since  the  elements  in  the  mass  matrix  are  constant,  { X ( s ) } 

-2  -2 
must  be  order  s  .  Therefore,  { X( s ) }  is  order  s  for  s  large  in  the 

right  hal f  s  plane. 

Having  now  established  that  response  of  the  structure  to  impulsive 

/\ 

loading,  {x(t)j, 

(x(t)}  =  L* 1  {X  (s)  }  (169) 
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exists  and  is  real,  continuous,  and  causal,  the  next  issue  to  be 
addressed  is  the  calculation  of  the  inverse  transform,  the  topic  of 
the  next  sectior. 
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SECTION  X 

CALCULATING  THE  RESPONSE  TO  IMPULSIVE  LOADING 


The  final  step  in  determining  the  impulse  response  of  the  structure 
is  to  calculate  the  inverse  transform  of  the  Laplace  transform  of  the 
response  to  impulsive  loading. 


(X(t)}  =  L  1  {X  (s)  } 


l-i  ?  un)Tu.)D(s) 

n-1  - 17 - TT“  {V 


mn(s  m-Xn^ 


The  inverse  transform  integral 


L1 {X  (s) } 


i  y+i“  t  . 

—  /  est  {X  (s ) }  ds 


z,wl  Y'10 


is  evaluated  using  the  residue  theorem  from  the  calculus  of  a  complex 
variable. 

The  closed  contour  of  integration,  used  in  conjunction  with  the 
residue  theorem,  is  given  in  Figure  18.  The  contour  is  divided  into 
six  segments  and  the  direction  of  integration  along  each  segment  is 
indicated  by  the  arrows.  Segments  3,  4,  and  5  of  the  contour  are 

1  /  m 

required,  because  the  branch  cut  and  branch  point  of  the  function  s 
are  taken  to  be  along  the  negative  real  axis  and  at  the  origin  of  the  s 
plane,  respectively. 

The  residue  theorem  states  that  the  integral  along  the  closed  contour, 
divided  by  2iri ,  is  equal  to  the  sum  of  the  residues  of  the  poles  of  the 
poles  of  the  integrand.  In  this  case,  the  statement  of  the  residue 
theorem  translates  into 

j-  /  (X(s)}es^  ds  =  -  I  I  ds  +  (172) 
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where  the  circled  index  indicates  the  contour  over  which  the  integration 
is  to  be  performed,  and  fcj  are  the  residues  of  the  poles  of  (X(s)} 
enclosed  by  the  contour. 


The  integrals  on  the  segments  of  the  closed  contour  are  evaluated 
for  the  case  where  the  length  of  segment  1  is  extended  indefinitely  in 
the  positive  and  negative  imaginary  directions. 

a  .  y + i  co  a 

/  {X(s)}est  ds  =  /  (X(s)}eSt  ds  (173) 

1  y-ico 

and,  as  a  result,  Equation  171  becomes 


1 

TH T 


Y+i"„ 

/  {X (s )  }e 


st 


ds 


y-i“ 


yiy  l  /  (X(s)}eStds  +  £b.  (174) 
k=2  i  J 


showing  that  one  need  evaluate  the  right  side  of  Equation  174  to  obtain 
the  inverse  transform. 


To  maintain  the  continuity  of  the  closed  contour,  the  radii  of 
segments  2  and  6  are  increased  indefinitely,  and  segments  3  and  5  are 
extended  indefinitely  in  the  negative  real  direction.  When  the  radii 
of  contours  2  and  6  are  increased  indefinitely,  it  can  be  demonstrated 
that  the  resulting  value  of  the  integrals  along  these  two  segments 
is  zero.  This  follows  directly  from  the  fact  that,  for  large  s,  (X(s)} 
is  order  s  .  Hence,  the  component  of  { X ( s ) }  is  order  s' 

+  -j  |s|  »  1  (175) 


Therefore,  the  asymptotic  expression  the  integral  on  segment  2  as  its 
radius  increases  is 


/  Xjl(s)estds 
2 


% 


lit 


Re" 

1 

Re  ° 


ds 


(176) 
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Expressing  s  in  polar  notation,  the  relation  is 


/  Xi(s)estds  *  j  ct  V-nrg-  Rie^de 
2  6o  R  e 


eRel8t  ie 


077) 


I s |  >>  1 


The  magnitude  of  an  integral  is  less  than  or  equal  to  the  maximum  value 
of  the  magnitude  of  the  integrands,  M,  times  the  length  of  the  path  of 
integration,  L. 


IT 

/  c 
9o 


eReiet 

*  1*7^ 


Reside 


<  M 


(178) 


The  maximum  value  of  the  magnitude  of  the  integral  is 

,Yt 


M  = 


ic*le 

R 


and  the  length  of  the  path  of  integration  in  radians  is 

L  =  n  -  eo 

The  resulting  bound  for  the  magnitude  of  the  integral  is 


M  •  L  = 


[111 

R 


(179) 


(180) 


(181) 


which  is  clearly  zero  for  finite  time  in  the  limit  as  R  becomes  large  or 
equivalently,  in  the  limit  as  s  becomes  large.  Since  the  magnitude  of  the 
integral  is  zero,  both  the  real  and  imaginary  parts  of  the  integral  are 
zero. 

lim  /  X(s)eStds  =  0.  +  iO.  (182) 

R-*-“  T 
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The  expression  for  the  magnitude  of  the  integral  along  segment  6  is 


0o  Re10t  ■ 

J  c  -  Re10ide 

-  R2ei2e 


<  M 


(183) 


and  the  proof  that 


lim  /  X  (s)estds  =  0.  +  iO , 

R--  6 


(184) 


follows  the  steps  given  in  Equations  179  through  181 


The  integral  along  segment  4  of  the  contour  is  evaluated  for  the  case 
where  the  radius  of  the  contour  shrinks  to  zero.  The  asymptotic 

t"  h 

expression  for  the  integrand,  the  £  n  component  of  Equation  163,  for 
small  s  is 


N  (s 

x,(s)  ^  l  -2 - _ 

i-1  mn(-»n) 


rni 


I s  |  <<  1 

111  «  u, 


(185) 


where  <p  p  is  the  £th  component  of  },  and 


D(s) 


1 . 


and 


I  s  1  <<  1 


(186) 


1 


(187) 
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The  asymptotic  expression  for  the  integral  on  segment  4  is 


/  X  (s)estds  -v  /  l 


7  — - r - +  estds 

r  7  ->  nl 

4  n=1  mn(-Xn) 


Is  I  <<  1 


is  <<  X. 


in  polar  notation,  where 


s  =  pe 


the  expression  the  integral  on  segment  4  is 


I  X  (s)e  ds  *  /  m(_-x  ) 

„  7i  n 

4 


N  {*  . _iet 

l  — - — 7 -  ep  tipe  6de  ( 1 90) 


Again,  the  magnitude  of  the  integral  on  segment  4  is  less  than  or  equal 
to  the  maximum  magnitude  of  the  integrand,  M,  times  the  length  of  the 
path  of  integration,  L. 

-n  N  U  }T{1. )  ie  - 
j  l  — 0 -  epe  tipe  6de  <  M  •  L  (191 ) 

*  n=1  mn('XJ 


The  maximum  magnitude  of  the  integrand  is 


N  U  )  U. } 

M  -  l  —2- - 

n=1 


e^l.p 


and  the  length  of  the  path  of  integration  is 


L  *  2  7T 
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The  resulting  expression  for  the  upper  bound  on  the  magnitude  of  the 
integral  on  segment  4  is 


/  X.(s)estds 


y  (*n}  tK} 
n=i  mn(-Xn) 


ept2ftp 


Is  |  «  1 

M  <<  Unl  (194) 

and  clearly  as  p  goes  to  zero,  or  equivalently  as  s  goes  to  zero,  the 
upper  bound  on  the  magnitude  of  the  integral  goes  to  zero.  Therefore, 

lim  /  X  (s)estds  =  0.  +  iO.  (195) 

P^O  4 

It  has  been  demonstrated  that  the  integrals  on  segments  2,  4,  and  6 
are  zero,  which  reduces  the  expression  for  the  inverse  transform  given 
in  Equation  174  to 

y+j.00 

—  J  {X  (s )  }eStds  =  --L.  I  /  {X(s)}eStds  +  £b.  (196) 

2iri  Y-i<»  2 tt i  k=3 , 5  ^  j  •* 

The  inverse  transform  is  seen  to  be  the  integrals  along  the  branch  cuts, 
segments  3  and  5,  plus  the  sum  of  the  residues  of  the  poles,  b.. 

J 

The  expression  for  the  integral  along  segment  3  is 

J  Xj,(s)estds  =  /PX£(reil,)erelirV1,dr  (197) 

3  R 


where 


s 


re 


i  it 


(198) 
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This  expression  simplifies  to 

/  Xt(s)estds 
3 


/  X4(re11T)e'rtdr 


(199) 


In  the  limit  as  the  radius  of  contour  2,  R,  goes  to  infinity  and  as  the 
radius  of  contour  4,  p,  goes  to  zero,  the  integral  along  segment  3 
becomes 

1  im  °° 

p-o  /  x?(s)estds  =  /  X  (re1 n)e *rtdr  (200) 

3  0  * 


Similarly,  the  expression  for  the  integral  along  segment  5  in  the  limit 
as  the  radius  of  contour  6  goes  to  infinity  and  the  radius  of  contour  4 
goes  to  zero  is 


lim 


p-*-0  J  X  (s)eStds 
R~  5 


X^  (re*ilT)e’rtdr  (201) 


The  sum  of  the  integrals  along  segment  3  and  segment  5  is 

.  I  I  X  (s)estds  =  /  «  Cre1,r)  -  X  Oe~iir))e'rtdr  (202) 

k  =  3 , 5  ,  *  0  *■ 

k 


Noting  that  Xf(re1TT)  and  X^(re~1TT)  are  complex  conjugates  of  each  other, 
the  expression  simplifies  to 


l  I  \(s)cstds 
k=3 ,  S  k 


2i  Im/  X  .  (re'11,)e'rtdr  (203) 


The  only  terms  in  the  inverse  transform,  (Equation  195),  remaining  to 

/*•  rf 

be  evaluated  are  the  residues  of  the  poles  of  X,,(s)e  .  By  definition, 

A  cf  ^ 

the  residues  of  the  poles  of  X^(s)e  are 


b. 

J 


1  im 
rm 

s^j 


,  Tm-iv  ,  St 

(S  -  Xj)X£(s)e 


(204) 
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or 


b. 

J 


lim  N  {*  }  {l.}D(s) 

s-xm  l  — ~ — 17 - 

j  J  "M  ™„(S  ) 

n  n 


<t>  e 
n 


st 


(205) 


which  reduces  to 


lim  (s-A™) 

2_ 


.  m 
s  -  A  . 
3 


(S  m-A,) 


U- K{l.}D(s)  t 

— - -  *  .  e 

m  3  l 

m  •  J 


(206) 


Before  taking  the  limit  to  evaluate  the  residue,  it  should  be  pointed 
out  that 


cs-;») 

l/m  r 
(s  -Aj 


4) 


m  -v  1  1  r/ 

1  s  m 

r=l  J 


(207) 


which  simplifies  the  expression  for  the  residue  to 


b. 

3 


lim 


m  ,  r/ 

J  a  T  ' 1  s 1  ‘  m 
:m  l  S  i  b 
s-*A j  \r-l  J 


(♦i )*(1. }D(s)  t 

— * -  4> .  e 

m  3  Jt 

m .  J 

3 


(208) 


Evaluating  the  limit  produces 


b.  . 


m . 
3 


>3*e 


Ajt 


(209) 


Having  evaluated  the  residues  of  the  poles  of  the  integrand  in 
the  inversion  integral,  the  evaluation  of  the  inverse  transform 
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is  complete.  The  general  expression  for  the  displacement  response  of 
the  structure  for  simultaneous,  unit,  impulsive  loading  at  all  its 
degrees  of  freedom  is 


x.(t)  =  — Im  [/  X  (re'l7T)e"rtdr] 

*  0 


mX1?'1!*.  }T{1.  IDCx1")  x“t 

+  l  - 1 - i-  *  e  J  (210) 

j  m.  J 


This  expression  is  based  on  Equation  195,  using  Equations  181,  203,  and 
209  for  the  appropriate  substitutions. 

Note  that  the  response  of  the  structure  has  two  parts,  one  part  being 
a  sum  of  decaying  sinusoids  and  the  other  part  an  integral  that  decreases 
with  increasing  time.  The  integral  does  not  decrease  exponentially, 
because  it  is  asymptotic  to  t"a  for  large  t,  where  a  is  greater  than  one. 
Therefore,  the  integral  dominates  the  response  for  a  time  long  after  the 
loading.  This  component  of  the  response  describes  the  non-oscillatory 
return  of  the  structure  to  its  unloaded  equilibrium  position. 

In  summary,  the  response  of  the  structure  to  impulsive  loading  can 
be  calculated  using  contour  integration  to  evaluate  the  inverse  transform. 
Having  obtained  the  response  to  impulsive  loading  as  a  function  of  time 
makes  possible  determination  of  the  response  of  the  structure  to  general 
loading  conditions  using  convolution.  In  essence,  the  calculation  of  the 
response  of  the  structure  is  no  longer  tied  to  the  use  of  Laplace  or 
Fourier  transforms.  In  other  words,  the  response  of  the  structure  to 
loading  time  histories  for  which  transforms  do  not  exist  can  be 
calculated  using  the  response  to  impulsive  loading  and  the  convolution 
integral . 
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SECTION  XI 

SUMMARY  AND  CONCLUSIONS 

The  generalized  der.vative  of  fractional  order  is  a  mathematical 
operator  well-suited  for  describing  the  frequency-dependent  mechanical 
properties  of  viscoelastic  materials.  As  shown  in  Section  V  and 
Appendix  C,  the  generalized  derivative  constitutive  relations  for  the 
materials  are  capable  of  describing  the  frequency-dependent  stiffness 
and  damping  of  the  materials  over  several  decades  of  frequency,  as 
observed  under  conditions  of  sinusoidal  motion. 

Moreover,  as  demonstrated  in  Section  V,  the  constitutive  relation 
for  3M-467  predicts  accurately  the  non-periodic  (transient)  response  of 
the  material  as  observed  in  the  laboratory.  In  addition,  the  generalized 
derivative  model  for  3M-467  performed  remarkably  better  than  the  Voigt 
model . 

The  generalized  derivative  RT  and  RTG  voscoelastic  models  enabled 
the  formulation  of  the  viscostiffness  matrix  of  the  viscoelastic  finite 
element.  This  led  to  the  successful  formation  of  the  equations  of  motion, 
for  a  structure  containing  both  elastic  and  viscoelastic  components, 
which  can  be  decoupled  and  solutions  obtained  using  modal  analysis. 

As  demonstrated  in  Section  VIII,  the  parameters  of  the  solutions  can  be 
determined  using  the  iterative  numerical  schemes  presented. 

Finally,  it  was  demonstrated  that  the  response  of  the  structure  to 
impulsive  loading  always  exists  and  is  continuous,  real,  and  causal. 

The  general  form  of  the  response  to  impulsive  loading  was  evaluated 
using  contour  integration. 

In  conclusion,  the  approach  to  viscoelasticity  resulting  from  this 
investigation  is  particularly  powerful,  in  that  the  general  motion  of  a 
structure  having  both  elastic  and  viscoelastic  components  can  be  determined 
given  that  two  conditions  are  met.  First,  the  RT  and/or  RTG  models  for 
the  viscoelastic  materials  in  the  structure  must  exist  and  comply  with  the 
second  law  of  thermodynamics  as  indicated  in  Appendix  A.  Second,  every 
component  of  the  structure  must  have  a  suitable  finite  element  approximation 
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APPENDIX  A 

THE  RT  AND  RTG  MODELS  AND  THE  SECOND  LAW 
OF  THERMODYNAMICS 

The  focus  of  the  thermodynamic  considerations  centers  on  the  ability 
of  the  RT  and  RTG  models  of  viscoelastic  materials  to  satisfy  the 
second  law  of  thermodynamics.  In  particular,  the  objective  is  to  derive 
constraints  on  the  parameters  of  the  RT  and  RTG  models  to  ensure  that 
the  state  of  stress  and  strain  in  the  material  predicted  by  the  models 
is  consistent  with  the  second  law.  The  considerations  are  limited  to 
the  case  where  the  material  is  at  a  constant,  uniform  temperature  and 
the  deformations  and  resulting  stresses  in  the  material  are  sinusoidal. 

The  assumption  of  constant,  uniform  temperature  is  motivated  by  the 
observation  that  the  moduli  of  viscoelastic  materials  are  usually 
strongly  dependent  on  temperature.  Allowing  the  temperature  of  the 
material  to  vary,  due  to  the  energy  dissipated  during  the  deformation 
process,  implies  that  the  parameters  of  the  RT  and/or  RTG  model  should 
be  modified  to  account  for  the  resulting  changes  in  material  properties. 
However,  the  RT  and  RTG  models  are  employed  where  the  parameters  of  the 
models  are  not  changed  during  the  deformation  process,  which  is  in 
essence  a  statement  that  the  temperature  of  material  does  not  change 
during  deformation. 

To  implement  the  assumption  of  constant,  uniform  temperature,  energy 
sinks  are  assumed  to  be  present  at  every  point  in  the  material.  The 
energy  sinks  instantaneously  absorb  all  the  dissipatfi  energy,  preventing 
the  temperature  from  changing. 

The  energy  sinks  are  represented  in  the  first  law  of  thermodynamics 

pe  =  p  +  V*FT  +  pq  (A.  1 ) 

by  q,  the  local  rate  of  change  of  energy  density  through  non-mechanical 
effects.  In  this  local  form  of  the  first  law,  V-h  is  the  divergence  of 
the  energy  flux  vector,  p  is  the  local  rate  at  which  internal  work  is 
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done  by  mechanical  means,  e  is  the  internal  energy  per  unit  mass  of  the 
material,  and  p  is  the  mass  density. 

The  following  development,  based  solely  on  the  first  and  second  laws, 
is  an  adaptation  of  the  work  of  Coleman  and  Truesdell  on  the  thermodynamics 
of  deformations  for  the  specific  case  at  hand  (Reference  29).  In 
notation  identical  to  Truesdell 's,  the  second  law  for  a  material  at 
constant,  uniform  temperature  is 

p6n  *  V*h  ■  pq  >  0  (A. 2) 

where  n  is  the  entropy  per  unit  mass  and  6  is  temperature.  This  local 
form  of  the  second  law  is  based  on  the  Clausius-Planck  inequality. 

Two  other  entities  relevant  to  the  following  discussion  are  the  free 

energy  of  the  material,  ip, 

ip  =  e  -  n6  (A. 3) 

and  the  internal  dissipation  rate,  6, 

6  =  p  -  p(n6  -  <0  (A. 4) 

where  the  dots  indicate  first  time  derivatives.  Substituting  the 
expression  for  the  dissipation  rate  produces 

6  -  p9n  +  (p  -  pe)  (A. 5) 

Using  the  first  law  to  substitute  for  p  -  pc  in  Equation  A. 5  yields 

6  =  pen  *  -  pq  (A. 6) 

From  the  second  law  and  the  above  expression  for  the  dissipation  rate, 
it  is  clear  that  the  second  law  is  always  satisfied  when 

6  >  0  (A. 7) 

In  other  words,  when  the  local  internal  dissipation  rate  of  a  material  is 
non-negative,  the  behavior  of  the  material  is  consistent  with  the  second 
law.* 

*This  result  is  contained  in  Coleman's  general  theorem  on  the  non¬ 
equilibrium  thermodynamics  of  deformation  (References  30,  32). 
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The  task  remaining  is  to  express  the  dissipation  rate  in  terms  of 
the  stresses  and  strains  in  the  material.  To  that  end,  recall  that  the 
energy  sinks  absorb  all  the  dissipated  energy. 

6  =  -pq  (A. 8) 

In  light  of  the  expression  for  the  dissipation  rate  given  in  Equation  A. 6, 

-V-fT  +  pen  =  0  (A. 9) 

Since  the  absorption  of  the  energy  by  the  sinks  prevents  the  opportunity 
for  energy  conduction, 

K  =  (A. 10) 


and,  as  a  result 


7*h  =  0 


and  consequently 


n  =  0 . 


(A. 11) 
(A. 12) 


Since  the  rate  of  change  of  the  entropy,  n,  is  zero.  Equation  A. 5 
reduces  to 

6  =  (p  -  pe) 


(A. 13) 


The  internal  dissipation  rate  can  now  be  expressed  in  terms  of  stress 
and  strain,  provided  the  internal  energy  per  unit  volume,  pc,  can  be 
expressed  in  terms  of  the  rate  of  change  of  strain  energy,  u.  The 
internal  energy  per  unit  mass,  c,  is  taken  to  be  a  function  of  temperature 
plus  the  strain  energy  per  unit  mass,  u/,p. 

e  =  f(e)  +  U/P  (A. 14) 

For  small  strain  the  expression  for  one  over  density  is 

1/p(t)  =  (1  +  e(t))/po  (A. 15) 

where  pQ  is  the  unstrained  density  and  e(t)  is  the  dilatation  strain, 
defined  in  Equation  38.  The  resulting  expression  for  the  internal  energy 
is 

e  «  f(e)  +  u*(l  +  e(t))/pQ  (A. 16) 
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which  simplifies  to 

e  =  f(0)  +  u/pQ 


(A. 17) 


for  small  strain.  Since  the  temperature  is  constant,  the  expression  for 
the  strain  energy  rate  is 


(A. 18) 


and  the  resulting  expression  for  the  dissipation  rate  is 

6  *  (p  -  u)  (A. 19) 

again,  for  the  small  strain. 


The  dissipation  rate  in  terms  of  the  stresses  and  strains  is 

6  =  ({ o  (t )  }T{  e  (t)  }  -  {o'(t)}T{e(t)})  (A. 20) 


or 

6  =  {  { o  (t )  }  -  (a ' (t)  }  }T{  e  (t)  }  (A. 21) 


where  f e( t ) }  is  the  column  vector  of  the  strain  rates,  (o(t)}  is  the 
column  vector  of  the  stresses  and  (o'(t)}  is  the  column  vector  of  those 
components  of  the  stress  acting  to  store  and  retrieve  strain  energy  in 
the  material.  As  indicated  in  Section  IV,  the  stresses  storing  and 
retrieving  strain  energy  for  steady-state,  sinusoidal  motion  were  those 
stresses  in-phase  with  the  strains.  Consequently,  the  total  stresses, 
{a(t)},  minus  the  in-phase  stresses,  (o'(t)},  leaves  the  out-of-phase 
stresses,  (o"(t)}. 

6  =  {o"Ct)}T{e(t)  }  (A. 22) 


The  stress  out-of-phase  with  the  strains 


e  „(t)  =  e  Sinu  t 

mn v  mnQ  o 


(64) 
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are  those  stresses  proportional  to  coswQt  in  Equation  65 

°"mn(t)  *  8mnF3(“o)eocos“ot  *  F4 t“o)t»T.0cos“ot  <A-23*> 

The  dissipation  rate  in  terms  of  the  stresses  and  strains  given  above  is 

3  3 

5  =  i»  F,(tu  )e2cos2w  t  +  l  I  a)  F.(w  )e2  cos2u  t  (A. 24) 

o  3  o  o  o  m_i  o  4^  o  mn„  o  '  ' 

n- l  m-l  o 

If  the  parameters  of  the  RT  and/or  RTG  models  are  chosen  such  that 


V"o>  5  0 

(A. 25) 

P,(«0)  2  0 

(A. 26) 

for  all  positive  w  ,  the  state  of  stress  and  strain  in  the  material 
predicted  by  the  models  is  consistent  with  the  second  law. 


Although  the  above  development  is  for  sinusoidal  motion  of  the 
material,  its  application  extends  to  periodic  strain  histories.  Consider 
the  case  where  the  material  undergoes  prescribed,  piece-wise  continuous 
and  bounded  strain  history  of  finite  duration  and  the  material  is  allowed 
sufficient  time  for  the  stresses  to  relax  to  their  unloaded  equilibrium 
value  of  zero.  Assume  that  at  some  later  time  the  same  strain  history 
is  imposed  on  the  material  and  the  material  again  allowed  to  fully  relax, 
and  this  process  of  loading  and  relaxation  is  repeated  continuously  at  a 
spc:ified  interval.  The  stress  and  strain  histories  of  the  material 
are  periodic  and  can  be  broken  into  their  respective  frequency  components 
using  Fourier  analysis.  The  internal  dissipation  of  each  frequency 
component  of  the  strain  is  given  by  Equation  A. 24  and  satisfies  the 
second  law,  given  that  Equations  A. 25  and  A. 26  hold. 


* F 3 ( <J,'0 ^  an<*  F4^wo^  are  functions  the  parameters  of  the  RT  and  RTG 
models  and  are  defined  by  Equations  67  and  70. 
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Consequently,  the  periodic  stress  and  strain  histories  of  the  material 
comply  with  the  second  law.  As  a  result,  satisfying  Equations  A. 25  and 
A. 26  is  sufficient  to  ensure  that  the  response  of  the  material  to  a 
piece-wise  continuous  and  bounded  strain  history  of  finite  duration 
satisfies  the  second  law.  This  follows  from  the  observation  that  each 
loading  and  relaxation  cycle  of  the  material  can  be  viewed  as  an 
independent  response  sequence. 


94 


AFML-TR- 79-41 03 


APPENDIX  B 

A  GENERALIZED  DERIVATIVE  RELATION  FOR  A 
NEWTONIAN  FLUID 


By  way  of  motivation  toward  describing  the  motion  of  physical  systems 
with  generalized  derivative  equations,  it  is  interesting  to  note  that  a 
Newtonian  fluid,  defined  by  the  stress-strain  rate  constitutive  relation 


o 


(B.l) 


where  u  is  the  bulk  viscosity,  can  undergo  dynamic  motion  which  lends 
itself  to  description  using  a  generalized  derivative.  When  the  Newtonian 
constitutive  relation  is  introduced  into  the  one-dimensional  momentum 
equation  for  a  homogeneous,  incompressible  fluid  at  uniform,  constant 
temperature 


3v  _  3o 
p  3t  3  z 


(B .  2) 


the  resulting  differential  equation  is 


p 


3  2v 
3Z2 


(B.3) 


This  equation  is  immediately  recognized  as  being  in  the  form  of  the 
one-dimensional  diffusion  equation.  Donaldson  (Reference  31)  has 
demonstrated  that  solutions  of  the  one-dimensional  diffusion  equation 
may  be  represented  using  fractional  calculus.  Of  particular  interest 
at  this  point  is  the  solution  to  the  one-dimensional  diffusion  equation 
that  occurs  when  an  infinite  half-space  of  such  a  Newtonian  fluid  is 
bounded  by  a  "wetted"  planar  surface  undergoing  some  known  motion.* 

If  the  known  motion  of  the  bounding  surface  has  a  Laplace  transform, 
Laplace  transforms  may  be  used  to  solve  the  diffusion  equation. 


*This  problem  is  a  generalization  of  Stokes'  second  problem  in  which 
Stokes  assumed  the  motion  of  the  "wetted"  plate  to  be  sinusoidal 
(Reference  33). 
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Figure  8-1.  Schematic  of  the  Half  Space  of  Newtonian  Fluid  Bounded 
by  a  "Wetted"  Surface 
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A  derivation  of  the  solution  is  given  here  for  completeness.  Taking 
the  Laplace  transform  of  Equation  B.3  results  in 

a  2V 

p(sV(s,z)  -  v  (0 ,  z) )  =  p  - —  (s,z)  (B.  4) 

az 2 

where  V(s,z)  is  the  transform  of  the  velocity 

oo 

V(s,z)  =  /  v(t,z)e  stdt  (B.5) 

0 

Assuming  the  fluid  is  initially  excited  from  a  state  of  rest,  the  initial 
velocity,  v(0,z),  is  zero  and  the  resulting  differential  equation  is 

3  2V 

psV  (s  ,  z )  =  -  (s  ,  z )  (B .  6) 

3  z  2 

The  solution  of  this  differential  equation  is  readily  seen  to  be 

/ P  ^  ^  / P  ^  ^ 

V(s,z)  =  A(s)e+  M  +  B(s)e  u  (B.7) 


Choosing  the  coordinate  system  such  that  the  normal  of  the  bounding 
surface  is  parallel  to  the  negative  Z  axis,  one  can  proceed  to  determine 
the  functions  A(s)  and  B(s)  in  terms  of  the  boundary  conditions.  The 
first  boundary  condition  to  be  satisfied  is  that  the  magnitude  of  the 
velocity  of  the  fluid  be  bounded  for  large,  negative  Z.  Consequently, 

B(s)  must  be  zero.  The  other  boundary  condition  is  that,  at  the  interface 
of  the  fluid  and  the  bounding  surface,  the  velocity  of  the  fluid  must  be 
equal  to  the  velocity  of  the  bounding  surface,  vQ(t).  To  enforce  this 
boundary  condition,  one  first  takes  the  transform  of  vQ ( t ) 

t 

Vn(5)  =  J  v  (t)e"stdt  (B . 8) 

o  o  o 


and  the  boundary  condition  is  applied  in  the  transform  domain  as  shown 


V  (s  ,  z ) 


(B.  9) 
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V(s  ,0)  =  A(s)  =  VQ(s)  (B.10) 


The  resulting  expression  for  the  transform  of  the  velocity  of  the  fluid 
at  a  distance  z  below  the  plate  is 


V(s  ,z) 


Vs)e 


(B.ll) 


At  this  point  one  wishes  to  determine  the  stresses  in  the  fluid 
according  to  the  Newtonian  constitutive  relation  (Equation  B.l).  The 
transform  of  Equation  B.l  is 


o*(s  ,z) 


3V  (s  ,  z ) 
3z 


(B.12) 


Substituting  Equation  B.ll  into  Equation  B.12  for  V(s,z)  produces 

o*(s,z)  =  /FF  /i“  V(s,z)  (B.14) 

The  transform  of  the  stress  may  be  inverted  by  observing  that  the 
transform  is  the  product  of  the  transforms  of  two  known  functions  of  time. 

o*(s,z)  =  7?  *  sV(s  ,z)  =  /ppL[(r(>j)tH)  X] 

•  L[v(t,z)]  (B.15) 

Thus,  the  stress  is  seen  to  be  the  convolution  of  the  two  functions  of 
time, 

o(t,z)  =  GL.  \  YJ 1*2*1  dt.  (B.16) 

r(>s)  o  (t-x)^ 

For  the  zero  initial  condition  on  velocity,  Equation  B.16  is  equivalent  to 

o(t,z)  =  A.  fZ  yJ±.>L\  dt  (B.17) 

r(»j)  dt  0  (t  ■  x )  * 
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Notice  that  Equation  B.17  states  that  the  stress  at  any  location  in 
the  fluid  is  equal  to  a  constant,  /up  ,  times  the  generalized  derivative 
of  fractional  order  ‘-2  of  the  velocity  of  the  fluid  at  that  location. 

o(t,z)  =  /u p  D*5  [v(t,z)]  (B.  18) 

This  stress-velocity  relation  evaluated  over  an  area  A  of  the  "wetted" 
surface,  z  =  0,  produces  a  force-velocity  relation. 

f(t,0)  =  A  /ST  Dh  [vo(t)]  (B .  1 9) 

Thus,  we  see  that  in  this  one  case,  the  macroscopic  behavior  of  a 
Newtonian,  viscous  fluid  is  characterized  by  a  generalized  derivative  of 
fractional  order  ‘-s,  even  though  the  microscopic  behavior  is 

a(t,z)  =  ue(t,z)  (B.  20) 

This  observation  suggests  that  generalized  derivatives  may  have 
applications  in  other  situations  where  a  global,  or  discrete,  description 
is  desired  of  a  phenomenon  which  is  locally  viscous. 
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APPENDIX  C 

RT  AND  RTG  MODELS  FOR  VISCOELASTIC  MATERIALS 

Generalized  derivative  constitutive  relations  for  three  viscoelastic 
materials  are  presented.  The  material  properties  on  which  the  RT  and  RTG 
models  are  based  {Reference  34)  were  determined  using  the  "temperature- 
frequency  superposition"  principle  (Reference  33). 

The  RTG  model  for  the  sheir  modulus  of  the  viscoelastic  material 
3M-467  at  75°F  is 

(1  +  b1D6l)o(t)  =  (u0  ♦  u1D“1)e(t)  (C.l) 

where 

bj  =  8  x  10~4  sec*  51  { C . 2) 

no  =  1.0  lb/in2  (C.3) 

=  7.3  lb-sec’ 56/in2  (C.4) 

=  .51  (C.  5) 

and 

Oj  =  .56  (C.6) 

The  mechanical  properties  predicted  by  the  model  are  compared  to  the 
material's  properties  in  Figure  C-l .  Note  the  excellent  agreement 
between  the  model  and  the  material  properties  over  8  decades  of  frequency. 


The  RT  model  for  the  viscoelastic  Young's  modulus  of  Sylgard  188  at 
1 20° F  is 


o(t) 


(Eq  ♦  E2D  1 ) e  (t) 


(C.7) 
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where 


E  =  60  lb/in2 

o 


Ej  =  43  lb-sec-  /in2 


(C.8) 


cij  =  .24 


(C.10) 


The  mechanical  properties  predicted  by  the  model  are  compared  to  the 
material's  properties  in  Figure  C -2. 

The  RT  model  for  the  viscoelastic  Young's  modulus  of  BTR  at  45°F  is 


o(t)  =  (eV1  +  E,D“2)e(t) 


(C.ll) 


where 


Ej  =  8S0  lb-sec*  /in2 


^2  =  18  lb-sec*  /in: 


~  *28 


(C.12) 


(2.13) 


(C.14) 


(C. 1 5) 


A  comparison  of  the  model  and  material  properties  is  given  in  Figure  C-3. 

Although  the  agreement  between  the  material  properties  and  their 
respective  models  is  very  good,  not  all  viscoelastic  materials  lend 
themselves  to  characterization  by  generalized  derivatives  models.  The 
materials  most  suited  to  modeling  with  generalized  derivative 
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constitutive  relations  are  those  that  have  properties  such  that  the 
following  relation  holds  in  the  transition  region: 

it  a  ,  . 

n  =  tan  —  (C.16) 

2 

where  n  is  the  loss  factor  and  a  is  the  slope  of  the  plot  of  the  log^ 
of  the  real  part  of  the  modulus  plotted  as  a  function  of  log-jg  of  the 
frequency  of  motion. 

In  summary,  generalized  derivative  constitutive  relations  do  in  fact 
model  the  frequency  dependent  mechanical  properties  of  at  least  three 
viscoelastic  materials.  However,  each  model  is  of  a  different  basic 
form,  as  can  be  seen  by  comparing  Equations  C.l,  C.7,  and  C.ll.  Hence, 
the  RT  or  RTG  models  are  capable  of  describing  viscoelastic  materials 
having  distinctly  different  mechanical  properties. 
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Figure  C-3.  The  Mechanical  Properties  of  BTR  Compared  to  the  RT  Model 
of  BTR 
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Loss  Factor  -  n(to) 
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